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SUMMARY 


A computer  simulation  of  the  signal  generated  by  SUS  has 
been  performed  for  the  purpose  o*  determining  the  accuracy  that 
could  be  achieved  in  the  performance  of  a source  level  experiment. 
It  is  postulated  that  for  1/3  octave  band3,  current  source  level 
estimates  are  correct  to  within  +2  dB  and  a new  experiment 
would  be  warranted  only  if  an  accuracy  of  +1  dB  could  be  achieved. 
With  the  exception  of  a few  problems  that  either  cannot  be  simu- 
lated, or  for  which  our  knowledge  of  the  physics  is  inadequate, 
it  is  concluded  that  an  accuracy  within  +1  dB  is  achievable. 

The  dif f iculities  arise  in  regard  to  the  yield  variability  of 
SUS,  the  form  of  non-linear itv  that  can  be  introduced  on  surface 
reflection,  and  the  transition  between  shock  arid  acoustic 
propagation  for  the  shock  wave  portion  of  the  signal.  Although 
the  latter  two  problems  have  been  included  in  the  simulation 
studies  on  the  basis  of  existing  knowledge,  a set  of  mini-experi- 
ments to  investigate  these  three  factors  is  recommended. 
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1.  INTRODUCTION 


SUS  shots  are  commonly  used  for  long  range  acoustic 

propagation  loss  measurements.  The  equivalent  acoustic 

source  level  of  explosives  has  been  known  approximately  for 

well  over  a decade"*".  As  acoustic  measurement  and  prediction 

accuracies  have  improved,  additional  measurements  and  analyses 

have  been  performed  by  many  scientists  to  more  accurately 

determine  these  levels.  An  advanced  prediction  model  in  common 

2 

usage  was  proposed  by  Gaspm  and  Schuler  in  1971.  Following 
the  CHURCH  ANCHOR  and  SQUARE  DEAL  Experiments  it  became  apparent 
that  different  segments  of  the  acoustic  community  were  still 
using  significantly  different  values  for  the  equivalent  acoustic 
source  level.  The  Manager,  LRAPP,  appointed  a committee  to 
study  this  problem.  The  committee  reported^  that  they  were 
unable  to  resolve  these  differences,  which  were  as  large  as 
six  decibels  in  some  frequency  bands,  and  recommended  the  per- 
formance of  a new,  carefully  controlled  experiment. 


D.  E.  Weston,  "Underwater  Explosions  as  Acoustic  Sources", 
Proceedings  of  the  Physical  Society,  76 , p 233-249,  1960. 

J.  B.  Gaspin  and  V.  K.  Schuler,  "Source  Levels  of  Shallow 
Underwater  Explosions",  Naval  Ordnance  Laboratory  Report, 
NOLTR  71-160,  October  1961. 

"SUS  Source  Level  Committee  Report",  XC  Report  112, 
prepared  for  LRAPP  by  Underwater  Systems,  Inc.,  November  1975. 


These  recommendations  were  the  subject  of  a SUS  Source 
Level  Workshop  held  in  July  1976,  and  reported  in  Ref  4.  The 
range  of  differences  in  SUS  source  levels  was  reduced  by  about 
3 dB  as  a result  of  newly  reported  information.  The  general 
consensus  of  the  workshop  was  that  SUS  source  levels  were 
known  to  about  + 2 dB  for  most  cases  of  interest  and  somewhat 
less  well  known  in  other  cases.  Further,  it  was  felt  that  by 
performing  a new  experiment  accuracies  within  + 1 dB  could  be 
achieved . 

> Subsequent  to  the  Workshop,  the  Manager,  LRAPP , established 

a working  group  to  recommend  a new  experiment  design  consistent 
. with  the  Workshop  conclusions.  Dr.  A.  F.  Wittenborn  served  as 

I chairman  of  this  working  group  and  Underwater  Systems,  Inc. 

| provided  technical  support.  Dr.  Wittenborn  has  requested  that 

Underwater  Systems,  Inc.  perform  a supportive  investigation  of  a 

i 

! 

i number  of  technical  problems  that  arise  in  the  performance  of 

measurements  and  the  processing  of  the  data.  In  essence,  the 
investigation  is  directed  toward  more  accurately  determining  the 

i 

accuracy  that  could  be  achieved  in  a new  experiment.  If  SUS 
source  levels  are  now  known  to  within  + 2 dB,  a new  experiment 
would  be  of  value  only  if  the  potential  accuracy  of  + 1 dB  esti- 
mated by  the  Workshop  could  be  achieved.  Additionally,  to  the 
extent  that  there  is  not  sufficient  current  knowledge  to  resolve 
. the  matter  for  all  factors  known  to  contribute  ,to  inaccuracy, 

Report  of  SUS  Source  Level  Workshop,  Airlie  Conference  Center, 
7-8  July  1976,  prepared  for  LRAPP  by  Underwater  Systems,  Inc. 

! and  Tracor,  Inc.,  July  1976. 
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Underwater  Systems,  Inc.  was  requested  to  recommend  a set  of 
mini-experiments,  each  of  which  would  be  designed  to  investigate 
a single  phenomenon.  Finally,  in  the  event  that  the  studies, 
and  mini-experiments  demonstrate  the  feasibility  of  achieving  an 
accuracy  of  +1  dB,  the  additional  knowledge  obtained  would 
permit  a reduction  in  the  size  of  the  experiment  that  would  have 
to  be  performed  to  determine  source  levels  within  that  accuracy. 

In  Chapter  2 of  this  report,  issues  that  had  previously 
been  raised  relative  to  the  sensor  and  electronic  systems  are 
addressed.  It  is  shown  that  a standard  steady-state  acoustic 
calibration  of  the  system  of  the  type  normally  performed  at 
Orlando  is  satif  actory.  Integration  of  the  signal  over  the 
sensor  extent  is  automatically  accounted  for  in  such  a calibration, 
provided  that  the  entire  system  is  linear. 

Chapter  3 is  concerned  with  signal  processing  and  data 
interpretation.  An  analytic  representation  of  an  explosive  signal 
has  been  developed  that  closely  replicates  the  description  given 
in  Ref  2 in  both  the  time  and  frequency  domain.  A closed  form 
Fourier  Transform  of  this  analytic  signal  is  obtained  that 
represents  "truth"  for  all  subsequent  operations.  The  effect  of 
the  measurement  range  on  equivalent  acoustic  source  level  is 
investigated  as  a guide  toward  determining  whether  it  is  possible 
to  perform  an  experiment  that  defines  the  range  beyond  which  the 
shock  wave  can  be  treated  as  an  acoustic  signal  within  tolerable 
error.  Digital  processing  is  examined  from  the  standpoint  of  the 
digitization  rate,  quantization  level,  registration,  and  system 
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bandwidth,  to  identify  the  errors  that  can  be  introduced  during 
signal  processing.  Finally,  the  practical  case  in  which  the 
received  signal  is  the  sum  of  the  directly-propagated  signal  and 
the  time-delayed  surface  reflection  is  simulated  for  both  linear 
and  non-linear  reflection  processes  to  investigate  procedures 
for  deconvolving  the  reflected  signal  and  determining  the 
errors  associated  with  this  process. 

Chapter  4 describes  a set  of  mini-experiments  to  investigate 
the  factors  for  which  our  basic  physical  knowledge  is  insufficient 
to  permit  a confident  estimate  of  the  errors  that  might  be 
introduced. 

Chapter  5 presents  the  conclusions  and  recommendations 


reached . 


2 . SYSTEM  RESPONSE 


The  two  major  factors  of  an  experimental  determination 
of  source  level  that  affect  the  overall  accuracy  of  the  exper- 
iment are  the  measurement  system  calibration  and  the  subsequent 
data  processing.  In  this  section  we  explore  three  factors 
which  impact  on  measurement  system  accuracy.  These  factors 
are: 

• Steady-State  and  Transient  System  Characteristics: 
as  chey  affect  system  calibration  and  calibration 
requirements. 

• Sensor  Size:  ef f .cts  on  the  system  calibration  and 
the  signal  that  is  acquired. 

• Calibration  Accuracy:  in  which  the  inherent  limita- 
tions of  various  calibration  techniques  are  stated. 

These  three  factors  are  examined  because  of  the  diverse 
views  held  by  the  members  of  the  various  committees.  It  is 
generally  agreed  that  small  sensors  and  electronic  systems 
which  have  a flat  response  with  zero  phase  shift  over  a wide 
.frequency  range  are  necessary  to  accurately  replicate  the 
pressure-time  history  of  the  signal;  e.g.  a small  sensor  is 
one  whose  physical  dimensions  are  small  relative  to  the  major 
features  of  the  explosive  signal.  It  is  further  agreed  that 
the  sensor-electronics  system  must  hive  an  adequate  dynamic 
range  and  a linear  pressure  response;  e.g.  the  system  must 


respond  linearly  to  both  the  high  levels  of  the  shock  wave 
and  the  much  lower  levels  of  the  negative  phases.  The  specific 
question  is  therefore  to  determine  the  acceptability  of  linear 
large  sensors  for  measurements  in  the  spectral  domain;  i.e.,  for 
use  in  determining  spectrum  levels. 

2.1  System  Steady-State  and  Transient  Characteristics 

The  frequency  and  phase  response  of  a data  acquisition 
system  is  known  to  "color”  (i.e.,  significantly  change)  the 
time  domain  waveform.  In  view  of  this  time  domain  coloration 
it  is  appropriate  to  ask  what  are  the  effects  in  the  frequency 
domain.  In  particular,  we  ask  how  frequency  and  phase  response, 
as  well  as  transient  response  of  the  measurement  system,  affect 
the  results  of  a measurement  of  the  energy  spectrum  of  a transient 
signal.  The  answer  to  this  question  is  widely  presumed  to 
be  that  the  squared  steady  state  frequency  response  is  the 
only  factor  that  affects  the  output  of  the  measured  energy 
spectrum,  but  we  find  no  straightforward  derivation  of  this 
fact  for  transient  signals  in  the  easily  accessible  literature. 

In  view  of  this  we  reproduce  here  such  a derivation. 

Consider  a time  signal  p(t)  which  has  zero  mean  and 
bounded  non-zero,  mean-squared  value  in  the  interval  0 to  T, 
and  is  identically  zero  outside  this  range,  i.e.. 
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p(t)  =0  for  t < 0,  T < t and 


and 


t+tc/2 


0 < 


J p2(t)  dt  £ M>  all  tQ  > 0,  0<^t^T 


t-tc/2 


Next  we  consider  applying  this  signal  to  the  input  of 
a linear  filter,  with  impulse  response  h(x),  frequency  response 
H (f ) (note  that  H(f)  is  complex).  h(x),  H(f)  are  a Fourier 


transform  pair  , i.e. 


H (u) 

h(x) 


J h(x)  e~jwTdx 
—00 

00 

— J H (w)  e^w^du) 
—00 


(2.1a) 


(2.1b) 


It  may  be  shown  that  for  a constant  harmonic  input  signal 
p^(iot)  = ae^wt,  the  output  signal  p0(u>t)  is  given  by: 

pQ(ujt)  = H ( uj ) p^  (wt) 

Thus  it  is  appropriate  to  call  H(w)  the  steady- state  frequency 
response,  since  it  describes  the  steady-state  change  in  phase 
and  amplitude  that  the  filter  applies  to  the  harmonic  input 
signal.  Our  input  signal  p^ft)  has  an  energy  spectrum  Ei(w) 
which  may  be  defined  from  the  time  signal  Fourier  transform, 
i.e, 




For  example  E.  B.  McGrab  and  D.  S.  Bloomquist,  The  Measure- 
ment of  Time  Varying  Phenomena,  pp  52-53. 
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Ei(a>)  = 2 Si(«)Si(«)  = 2|Si(w) 


where 


w 

S^w)  = f pi(t)  e~ja)tdt 


(2.2) 


(2.3) 


It  is  this  energy  spectrum  that  we  wish  to  measure.  The  out- 


put of  the  linear  filter  may  be  represented  as: 


00 

P0(t)  = f h(T)  s^t-x)  dT 


(2.4) 


or,  applying  the  Fourier  transform  (Eq.  2.3) 


SQ  (u)  » / / h(t)  sjL(t- 


t)  dx  e ^utdt 


-00  —00 


I h(Ti 


-3WT 


(2.4a) 


go 

J si(t-T)  e~ju>(t"x)dt  dx  (2.4b) 


J h(x)  e"^Tdx  S.(w) 


(2.4c) 


= H(<i))S^(uj) 


Thus ; 


E0(w)  « 2 j SQ  (u>) 


* * 


« 2(S.  <uj)  H (uj)  Si  (w)H  (u)  ) 


2 2 j j 

2 |s^  (u>)  j i H (u>)  I a E.(w)  |H(w)  I 


(2.5) 


We  note  that  the  frequency  response  function  can  be  sepa- 
rated into  magnitude  and  phase  terms,  i.e.: 


A (w)  ei9  ^ where  A(u»)  is  real. 


Then? 
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A2  (u>) 


E^(u))  = Eq(u))/A  (w)  (2.6) 

Thus  the  magnitude  of  the  steady-state  frequency  response  is 
a sufficient  calibration  of  the  (linear)  mo-^urement  system 
for  determining  the  input  energy  spectrum  at  a frequency  range 
of  interest#  if  the  output  energy  spectrum  from  the  system 
is  known. 

2.2  Sensor  Size  Effects 

As  in  the  case  of  using  steady-state  frequency  response 
to  characterize  a system's  response  to  a transient  signal, 
it  is  also  appropriate  to  ask  how  the  finite  size  of  a transducer 
affects  the  system's  ability  to  measure  a signal  whose  spatial 
extent  in  the  water  is  comparable  to  or  smaller  than  the  trans- 
ducer dimensions.  In  this  case,  we  will  show  for  a linear 
measurement  system  that  the  pure  tone  frequency  response  cali- 
bration (ratio  of  squared  voltage  to  squared  pure  tone  pressure 
as  a function  of  frequency)  in  the  direction  of  incidence  is 
a sufficient  calibration  to  permit  the  computation  of  the  exact 
energy  spectrum  for  the  incident  wave  form  in  the  frequency 
range  of  observation. 

Physically,  this  problem  is  sometimes  analyzed  as  an 
integration  of  the  sound  field  over  the  transducer  surface, 
so  that  the  transducer  output  for  an  arbitrary  input  has  the 


form? 
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p (t) 


/ P<*» 


t)  dS 


where  p(x,t)is  the  time  and  space  varying  pressure  field 

I 

p (t)  is  the  "apparent  pressure"  that  a point  trans- 
ducer would  see  in  order  to  produce  the  same 
output 

SQ  is  the  active  area  of  the  transducer 
dS  is  differential  element  of  surface  area 
A more  useful  representation,  which  can  be  used  when  the  sound 
field  has  a particular  direction  of  incidence  say  xQ,  is  the 
impulse  response  function,  defined  byj 


Vo(t)  = / h(T,XQ)  p(t-T,XQ)dT 


(2.7) 


The  impulse  response  function  h(T,xQ)  is  the  voltage  response 
of  the  transducer  to  a unit  amplitude  plane  wave  pressure 
impulse  arriving  from  the  direction  xQ.  This  formulation  is 
equivalent  to  the  former  representation,  where  the  integration 
over  the  transducer  surface  has  been  replaced  by  an  integration 
over  the  time  variation  of  the  (plane  wave)  pressure  wave  form. 
Taking  the  Fourier  transform  as  in  Eq.  2.4a  yields; 


Vo(o) 


J vo(t)e~^ajtdt 


= H(w,Xq)  S^(w,xq) 


(2.8) 


1C 


i 


where:  H(u),xq)  = 


j h(t,x0)e  J“'wdt 


is  again  the  complex  steady-state  frequency  response 
of  the  transducer  (i.e.  the  output  for  a constant 
input  sinusoid  of  unit  amplitude,  frequency  u) . 
Si(w,xQ)  is  the  complex  frequency  spectrum  of  the 
input  pressure  signal. 

The  energy  spectrum  of  the  output  voltage  is; 


Ev(w)  = 2 V0(w)  vo*<o>) 


= |H(w,X0)|  Ep  (to) 


where  [ H (u> # xq)  | = w(2Ttf,xo) 

is  the  transducer  sensitivity  (squared  voltage  response) 
as  a function  of  frequency  for  the  direction  xQ. 

Thus  we  see  that  the  steady-state  frequency  response  provides 
the  only  calibration  necessary  for  determining  energy  density 
of  an  incoming  plane  pressure  wave  of  arbitrary  time  dependence, 
Physically,  this  occurs  because  exactly  the  same  integration 
process  occurs  for  each  frequency  in  the  calibration  as  occurs 
for  the  respective  Fourier  frequency  component  in  the  incoming 
arbitrary  wave  form.  The  only  assumptions  necessary  to  prove 
this  were  that: 

1.  The  incoming  signal  be  preserved  for  the  entire 
time  during  which  it  is  non-zero. 

2.  The  transducer  is  linear  and  its  sensitivity  is 
known  as  a function  of  frequency  for  the  frequency 
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range  of  interest  and  the  directions  of  incidence 
of  the  acoustic  signal. 

3.  The  signal  consists  of  a plane  or  spherical  wave 
whose  direction  corresponds  to  that  for  which  the 
calibration  is  known. 

The  latter  requirement  concerning  the  spatial  form  of  the  wave 
front  is  necessary  in  order  to  uniquely  assign  an  impulse 
response  function  for  Eq.  2.7,  unless  the  transducer  is  omni- 
directional . 

2.3  Transducer  Calibration  Accuracy 

The  two  major  transducer  types  of  interest  are  tourmaline 
gauges  and  hydrophones. 

The  tourmaline  gauge  and  its  calibration  are  described  by 
Cole.®  He  cites  comparability  of  static  and  reciprocity  cali- 
brations of  within  5%  or  less  than  0.5  dB.  Bobber7  describes 
hydrophone  calibration  procedures,  of  which  the  reciprocity 
method  is  considered  the  most  accurate.  Reproducibility  of  these 
calibrations  is  again  typically  within  0.5  dB. 

Of  potential  significance  is  the  fact  that  neither  calibration 
is  typically  performed  in  the  pressure  regime  of  interest  for 
this  study  (eg.,  hydrophones  are  typically  calibrated  at  signifi- 
cantly lower  levels,  while  tourmaline  gauges  are  usually  calibrated 

y — — 

R.  H.  Cole,  Underwater  Explosions,  Dover,  New  York,  1948 
Section  5.7 

7 R,  J.  Bobber,  Underwater  Eiectroacoustic  Measurements,  Naval 
Research  Laboratory,  Government  Printing  Office, 

Washington,  D.C.,  1970 
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at  much  higher  levels).  As  a result,  dependence  on  gauge  and 
system  linearity  will  be  required  in  order  to  obtain  the  desired 
data. 

In  either  case,  a calibration  accuracy  for  the  transducer 
of  less  than  0.5  dB  appears  to  be  within  reach. 


3.  A SUS  SIGNAL  PROCESSING  MODEL 

The  analysis  of  transient  signals  using  digital  proces- 
sing techniques  has  proceeded  on  the  basis  of  analogy  with 
the  analysis  of  continuous  signals.  We  are  not  aware  of  any 
systematic  study  of  errors  that  may  occur  in  the  analysis 
which  take  into  account  some  of  the  unique  properties  of  tran- 
sient signals  (such  as  thei*-  bounded  total  energy). 

In  this  section  we  present  a comprehensive  analysis  of 
signal  processing  error  as  a function  or  the  processing  param- 
eters so  as  to  place  quantitative  estimates  on  those  error 
sources.  Other  sources  of  variability  or  error  in  an  actual 
experiment , including  variation  of  the  medium  and  background 
noise  are  not  included. 

The  digital  signal  analysis  process  is  represented  as 
consisting  of  three  distinct  processes,  which  include: 

• Sampling  in  the  time  domain  (temporal  sampling) 

• Quantization  of  the  time  domain  signal  into 
discrete  levels  (quantization) 

• Discrete  Fourier  Transform  (OFT) 

Each  of  these  processes  may  introduce  errors  that  reflect 
themselves  as  differences  between  the  resulting  estimate  of 
the  signal  energy  spectrum  and  the  "true*  spectrum  that  would 
result  from  performing  a continuous  analytic  Fourier  transform 
on  the  continuous,  exact  signal.  We  presume  for  the  purposes 
of  this  study  that  these  errors  may  individually  be  made  suf- 
ficiently small  that  their  interactions  are  negligible,  and 
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Analytic  Model  of  Explosion 
(Parameters:  yield,  depth, 
range,  reflection  delay, 
reflection  propagation  range) 

(See  Appendix  A) 


Fourier  Transform  of 
Analytic  Model 

SHOTSPECTRUM  II 


Sampled  time  series  of 
analytic  anti-alias  filtered 
signal- (sample  rate,  cutoff  freq.) 

SHOTTIME  II 


Quantization  of  signal  amplitude 
(scale  factor) 


Discrete  Fourier  transform 


FFTANAL 


Examine  data  differences 
in  narrowband  and  1/3  octave 
formats 

DATADIF 


Figure  3-1. 


Output  comparison 
graphs 


Schematic  block  diagram  of  computer  experi 
caents,  identifying  component  programs. 


the  total  mean-square  error  is  the  sum  of  the  component  mean- 
square  errors. 

This  section  includes/  where  possible,  analytic  estimates 
of  an  upper  bound  to  the  error  for  specific  error  sources. 

These  error  estimates  are  then  compared  with  empirical  error 
estimates  derived  from  computer  based  experiments.  The  experi- 
ments are  based  on  comparing  the  output  of  a simulation  of 
a digital  signal  processing  system,  operating  on  a sampled 
analytic  signal,  with  an  exact  Fourier  analysis  of  that  same 
analytic  signal. 

Figure  3-1  presents  a schematic  diagram  of  the  data  flow 
in  the  computer  experiments  reported  in  this  section.  The 
concept  of  the  experiment  is  to  use  a superposition  of  ana- 
lytic functions  (exponentials,  half-sines  and  constants)  to 
represent  the  time  series  of  an  explosive  signal.  This  func- 
tion may  be  transformed  analytically  to  obtain  an  exact  expres- 
sion for  the  energy  spectrum  of  the  time  series.  This  exact 
transformed  analytic  expression  is  evaluated  at  discrete  fre- 
quencies in  the  program  SHOTSPECTRUM.  SHOTSPECTRUM  permits 
the  computation-of ' a spectrum  for  a single  direct  arrival, 
or  the  superposition  of  a direct  arrival  plus  a surface-reflected 
arrival  with  arbitrary  delay,  propagation  range  and  amplitude 
clipping  (of  the  reflected  shock  wave).  The  program  SHOTTIME 
piovides  the  time  series  values  of  the  analytic  model  signal 
after  it  has  been  operated  upon  by  an  analytic  representation 
of  a 4-pole  Butterworth  anti-aliasing  filter.  Input  data  to 


this  program  includes:  sampling  time  interval,  shot  yield, 
depth  and  range,  filter  cutoff  frequency,  reflected  signal 
delay,  range  and  shock  amplitude  clipping.  The  program  FFTANAL 
operates  on  the  output  of  SHOTTIME  to  quantize  the  amplitude 
representation  of  the  signal,  and  then  implements  a Fast 
Fourier  Transform  (FFT)  to  provide  the  digital  signal  proces- 
sor estimate  of  the  signal  spectrum.  Finally,  the  program 
DATADIF  takes  the  output  from  SHOTSPECTRUM  and  FFTANAL  and 
plots  the  respective  spectra  and  their  differences  in  narrow- 
band  and  1/3  octave  formats  for  the  frequency  ranges  of  0 - 
100  Hz  and  0 - 500  Hz. 

There  are  a number  of  different  questions  which  can  be 
examined  using  this  model.  As  in  the  previous  chapter,  the 
questions  selected  for  investigation  are  those  which  could 
not  be  definitively  resolved  by  the  various  committees  and 
Workshop.  These  questions,  each  followed  by  a short  discus- 
sion are: 


• What  accuracy  can  be  achieved  with  FFT  processing? 

A comparison  of  spectral  levels  obtained  by  several  organ- 
izations by  applying  FFT  techniques  to  a set  of  SUS  record- 
ings is  reported  in  Ref.  3.  Differences  of  the  order 
of  a decibel  were  observed,  but  were  not  investigated 
in  detail,  since  the  goal  was  to  attempt  to  explain 
reported  differences  of  the  order  of  3 to  6 dB.  Since 
this  difference  has  since  been  reduced  by  3 dB,  Ref.  4, 
the  remainina  reported  processing  differences  become 
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important,  particularly  in  regard  to  the  performance  of 
an  experiment.  A series  of  investigations  are  therefore 
described  which  indicate  the  level  of  accuracy  which 
can  be  achieved. 

• What  accuracy  limitations  are  imposed  by  overlap 

of  direct  and  surface-reflected  signals?  Received  signals 
of  SUS  detonations  include  both,  direct  and  surface-reflected 
arrivals.  For  60  foot  1.8  lb  SUS  it  is  not  possible  to 
avoid  overlap  of  these  signal  components.  The  decorrela- 
tion of  the  received  signals  to  determine  the  free  field 
source  levels  is  examined  in  the  frequency  domain.  This 
investigation  includes  consideration  of  non-linear  effects 
at  the  surface. 

• What  is  the  optimum  measurement  range?  It  is  well 
known  that  the  shock  wave  generated  in  a detonation  propa- 
gates by  shock  laws  rather  than  acoustic  laws.  The  latter 
govern  the  other  parts  of  the  total  signal.  Consequently 
the  signal  spectrum  is  range  dependent.  It  is  believed 
that  as  the  propagation  range  increases,  and  the  shock 
wave  amplitude  decreases,  a gradual  transition  between 
shock  wave  and  acoustic  wave  propagation  occurs.  Further, 
it  is  believed  that  this  occurs  in  the  vicinity  of  a peak 
pressure  of  about  1 psi.  This  suggests  that  if  the  source 
levels  are  measured  at  sufficient  range,  they  will  ade- 
quately represent  the  "effective  source  level"  for  acous- 
tic measurements  at  greater  range,  to  within  a fraction 


of  a decibel.  However,  as  the  range  increases  the  accuracy 
of  refraction  corrections  to  propagation  decreases  and 
for  that  reason  one  would  want  to  make  the  source  level 
measurement  at  the  shortest  possible  range.  Accepting 
a peak  pressure  of  1 psi  as  being  approximately  correct 
for  the  transition,  an  investigation  of  potential  errors 
is  made  in  the  vicinity  of  this  value  to  determine  whether 
it  is  desirable  to  experimentally  account  for  our  limited 
knowledge  by  making  measurements  at  several  ranges. 

3.1  Analytic  Representation  of  the  Signal 

As  previously  discussed,  the  investigation  is  performed 
by  establishing  an  analytic  signal  which  closely  replicates 
the  SUS  signal  but  for  which  the  Fourier  Transform  can  be 
taken  in  closed  form.  No  attempt  has  been  made  to  develop 
an  analytic  model  which  represents  "truth"  in  the  absolute 
since.  Rather  it  is  sufficient  for  our  purpose  to  closely 
replicate  the  recognized  features  of  explosive  signals  in 
the  time  and  frequency  domains.  Thus,  the  analytic  pressure- 
time model  and  its  Fourier  Transform  are  used  to  represent 
a "relative  truth"  against  which  all  other  analyses  are  com- 
pared to  permit  the  systematic  development  of  error  estimates. 

Details  of  the  analytic  model  and  these  programs  are 
described  in  Appendix  A.  Figures  3-2  and  3-3  provide  plots 
of  the  time  series  signal  for  one  sat  of  shot  parameters, 
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showing  the  difference  between  the  wide  band  signal  and  the 
anti-alias  filtered  signal.  Figures  3-4,  3-5  and  3-6  show 
comparison  spectra  for  this  analytic  model  and  other  models. 
Of  principal  interest  here  is  not  detailed  agreement,  but 
rather  the  fact  that  the  time  series  used  in  this  study  pro- 
vide spectra  that  incorporate  many  of  the  features  of  the 
spectra  to  be  expected  from  other  models  or  from  an  actual 
experiment.  Thus  we  expect  that  error  analysis  estimates 
made  on  the  basis  of  studying  this  signal  should  be  closely 
representative  of  what  may  be  expected  during  analysis  of 
field  data. 


3.2  Digital  Processing 


3.2.1  The  Discrete  Fourier  Transform.  The  Discrete  Fourier 
Transform  (DFT)  is  the  analytic  construct  that  forms  the 
basis  of  modern  digitally- implemented  signal  processing.  It 
is  a close  relative  of  the  well  known  Integral  Fourier  Trans- 
form, but  has  some  very  distinctive  properties  o*  its  own. 

As  has  been  noted  previously  in  this  report,  much  of  the  present 
understanding  of  the  significance  of  the  DFT  is  based  upon 
the  study  of  spectral  analysis  of  continuous,  stationary,  in- 
finite  time  series.  In  the  application  to  noise  time  series, 
the  squared  magnitude  of  the  DFT  provides  an  estimate  of  the 
power  spectrum  of  the  infinite  length  time  series.  We  will 
show  that  the  DFT,  when  applied  to  the  entire  time  interval 


For  example  see,  "R.  B.  Blackman  and  J.  W.  Tukey,  The  Measure- 
ment of  Power  Spectra,  Dover,  New  York,  1968. 
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3-3.  Antialias  filtered  time  series.  (Cutoff  frequency  = 2040  Hz),  all 
parameters  same  as  Figure  3-2. 


Frequency,  Hz 

Figure  3-5  Comparison  of  USI  model  with  NSWC  model  - 1,8  lb  of  TNT  at 


Figure  3-6  1/3  octave  band  source  level  comparison  with  NSWC 

1.0  lb  of  TMT  at  60  ft  depth. 
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of  the  finite  duration  transient,  yields  the  exact  Fourier 
transform  of  the  signal  at  the  discrete  frequencies  at  which 
it  is  evaluated,  and  further,  when  sampled  at  sufficiently 
high  rate,  provides  sufficient  information  to  define  the 
continuous  Fourier  transform  at  all  values  of  frequency. 

Let  us  introduce  the  discrete  operator  (after  Blackman 

g 

and  Tukey,^  paragraph  A,2) 

00 

d (t)  = l At  6 (t-qAt) 
q=*-® 

which  has  the  Fourier  transform, 

00 

/»  oo 

d (t)  e“32lTft  dt  = l 6 (f-q/ At) 


where:  At  is  a non-zero  time  interval 

q is  an  integer  series  index 
6 the  Dirac  delta  function  defined  by; 


/ 


\f  (t)  A (t-c)dt  => 


f(c)  if  a < c < b 

0 if  c < a or  b < c 


Now,  if  we  take  the  disciete  Fourier  transform  of  a continuous 
bounded  signal  p(t)  we  have, 

an 

d(t)  p(t)  dt 


/ 


00  00 


l f At  p(t)  6 (t-qAt)  e“}2nft  dfc 

qO-®  J 

«rO> 

• l At  p(qAt)  e“32nf^t  » p'(f) 


q 


9— Co 
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We  name  this  function  P (f ) . From  the  exact  Fourier  transform 
of  p(t) , i.e. , 

OD 

P(f)  - J p(t)  e“j2lT?t  dt 

—a > 

l 

we  observe  that  P (f)  is  the  convolution  of  P(f)  with  D(f), 
(the  Fourier  transform  of  d(t)),  according  to  the  multiplica- 
tion theorem.  Thus, 

p'(f)  « P (f ) *D (f ) 


P(f-X)  D (X)  dX 


= l P(f  - q/At)  (3.2) 

qrs— a> 

Equation  3.2  is  a complex  but  extremely  significant  result. 

Let  us  first  examine  its  significance  for  the  case  where  the 
signal  has  no  energy  above  a frequency  fN  = l/2At  (the  Nyquist 
frequency)  i.e.; 

P(f)  » arbitrary,  f < fs 
°'  f i fN 

% 

Figure  3-7  is  a plot  of  the  resulting  real  part  of  P (f) 

resulting  from  the  transform  of  a real  time  series.  The  real 

♦ 

part  of  P(f)  is  an  even  function.  P (f)  is  a continuous  even 
function  of  frequency  f,  which  is  identical  with  P(f)  in  the 
region  < f < It  is  repeated  an  infinite  number  of  times 

centered  at  the  frequencies  +q/At  where  q » 1,  2 ... 
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figure  3-7  Comparison  of  discrete  and  integral  Fourier  transforms.  The  real  part 
of  a discrete  Fourier  transform  (Re(P'))  is  shown  plotted  below  the  real  part  of 
an  integral  Fourier  transform  (Re(P)>  for  a time  series  with  zero  spectral  enerqv 
components  at  frequencies  f greater  than  f„  = l/2&t 


Thus  our  condition  that  P(f)  have  zero  value  above  f„,  is 

N 

a sufficient  condition  for  P (f)  = P(f),|fj<f  . Further,  it 
is  a necessary  condition,  since  if  P{f)  / 0 for  some  f^  qreater 
than  fN  but  less  than  2f,  then  at  the  frequency  f2  = 2f^  - 
we  have  from  Equation  3.1, 


P (f2) 


= P(f2)  + P <f2 
» P{f2)  + P(f2 
= P(f2)  + P(fx) 
+ P<f2) 


- 1/At) 


- 2 f 


max 


) 


The  influence  of  signal  above  the  Nyquist  frequency  on 
the  resulting  DFT  is  called  aliasing,  and  its  removal  is 
called  anti-alias  filtering.  Thus  we  see  that  appropriate 
anti-alias  filtering  is  both  necessary  and  sufficient  to 
produce  a DFT  that  is  exactly  equal  to  an  I FT  of  a finite 
dt’-^tion  signal. 

When  we  perform  a Fast  Fourier  Transform  (FFT) , this 
amounts  to  an  evaluation  of  a DFT  (Equation  3.1)  at  discrete 
frequency  intervals  that  are  sub-multiples  of  f , if  the 

total  number  of  frequencies  is  N (usually  a power  of  2),  we 

* 

evaluate  the  complex  value  of  P (f)  at  the  frequency  series, 
n f 

f = n ° o*  l*  2f***N~i* 

9 

It  may  be  shown  that  a necessary  and  a sufficient  sampling 
for  completely  recovering  the  energy  spectrum  is  obtained  when 


^ w.  T.  Cochran  et  al,  "Burst  Measurements  in  the  Frequency 
Domain",  Proceedings  of  the  IEEE  Vol.  $4,  No.  6,  pp  830-841 
Appendix  9.2. 
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the  spectrum  is  computed  at  the  frequency  intervals  1/2T,  where 
T is  the  total  length  of  the  time  series.  For  a burst  dura- 
tion 0.5  seconds  long,  this  requires  spectrum  sampling  at 

every  1 Hz.  If  such  a sampling  is  achieved,  according  to  the 

9 3 

Cardinal  Series  of  Interpolation  Theory  the  continuous  Fourier 
spectrum  can  be  constructed  from  the  expression; 


P (f) 


l P’(n/t) 
n=~°° 


sin(7r(f  - n/T) T) 
7T(f  - n/T)  T 


(3.3) 


where  f is  an  arbitrary  frequency  in  the  range 
n is  an  integer  series 
T is  the  total  period  of  the  signal 


f 


<f 


max 


Note  in  the  above  expression  that  although  P repeats  cycli- 
cally at  frequencies  2nf  _ . the  sin  x/x  term  rapidly  damps 
out  so  that  the  computation  can  be  achieved  to  an  arbitrary 
level  of  accuracy  with  a finite  number  of  terms. 

Further,  the  energy  density  can  be  explicitly  evaluated 
9b 

also,  using  a similar  interpolation  expression,  i.e.. 


E(f)  - | P (f ) 

I 1 2 

where  |p(n/2T)| 


1 2 ® , 

I = l j P (n/2T) 

n=-«o 

are  discrete  samples  of  the  contiguous  energy 


sin  ( 2 fr  (f  - n/2T) T) 
2V(f  - n/'llT't 


(3.4) 


density  spectrum  jP(f)|  . 

Thus  we  see  that  the  DFT  is  an  extremely  powerful  way 
of  evaluating  impulse  spectra  with  the  possibility  of  allowing 


9 ci 

Cochran,  op.cit,.  Appendix  B.l 
9b 

Cochran,  op.cit. , Appendix  D 
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us,  in  a single  evaluation  of  the  data,  to  determine  exactly , 
all  of  the  information  that  is  available  in  the  energy  spec- 
trum, which  can  then  be  used  to  compute  exactly  the  energy 
spectrum  in  any  other  bandwidth.  The  requirements  for  this 
exact  analysis  using  any  DPT  are  the  following: 

• Th®  time  series  to  be  analyzed  contains  the  entire 
event. 

• The  energy  spectrum  is  computed  at  frequency  inter- 
vals given  by  Af  = 1/2T,  where  T is  the  duration  of  the 
event.  This  can  be  achieved  by  adding  zeros  so  that 
the  total  signal  length  analyzed  is  at  least  2T.  (This 
permits  application  of  the  interpolation  formula.) 

• The  energy  spectrum  of  the  time  series  contains  no 
energy  at  or  abo^e  a maximum  frequency  f = l/2At  (this 
elimates  all  aliasing  error) . 

The  third  condition  is  virtually  impossible  to  realize 
in  practice,  so  ir  the  following  section  we  derive  error 
estimates  for  realizable  anti-aliasing  filters,  after  a brief 
discussion  of  errors  arising  in  the  FFT  computation  due  to 
finite  word  length. 


3.2.2  Finite  Word  Length  FFT  Error.  The  procedure  almost 
universally  used  to  perform  a Discrete  Fourier  Transform  (DFT) 
is  the  "fast  fourier  transform"  (FFT) , which  is  a particular 
means  of  coding  the  DFT  operation  in  a way  that  substantially 
reduces  computation  time  over  that  of  direct  computation  by 
taking  advantage  of  periodicity  in  the  computation  procedure. 

It  is  inherently  no  more  nor  less  accurate  than  any  other 
means  of  computing  the  desired  result.  The  fact  that,  in  any 
computer,  the  numerical  values  used  must  be  of  finite  length 
implies  that  in  the  process  of  performing  the  many  additions 
and  multiplications  of  data  values,  the  rounding  or  truncation 
that  must  be  performed  at  various  levels  will  introduce  errors 
in  the  computed  results,  independent  of  the  errors  in  the  in- 
put data  to  the  program.  In  some  sense  then,  the  errors  may 
be  viewed  as  a noise  term  that  is  introduced  within  the  FFT 
program.  Unfortunately  the  magnitude  of  this  noise  is  far 
from  independent  of  the  input  signal.  Extensive  analyses  of 
FFT  error  have  been  published  in  the  literature,  and  will  not 
be  reproduced  here,  but  their  results  will  be  quoted.10 

For  a fixed  point  FFT,  where  the  data  values  are  presented 
as  integers  and  a special  form  of  rounding  is  used  (rounding  a 
value  up  or  down  at  the  edge  of  the  round-off  interval  is  per- 
formed on  a random  basis)  results  in  an  upper  bound  error  of, 

e.  R.  Rabiner  and  C.  M.  Rader,  ed.,  Digital  Signal  Processing, 
IEEE  Press,  N.Y.„  1972.  See  in  particular : P . D.  Welch,  “A 
Fixed-Point  Fast  Fourier  Transform  Error  Analysis",  pp  470-476, 
C.  L.  Weinstein,  "Round  Off  Noise  in  Floating  Point  Fast 
Fourier  Transform  Computation",  pp  477-483. 

33 


*!“•£ ;-?*•** 

jfcMric  ii  T i"  i i ■ £l iL £.  ^ i rii 


rms (error) 
rms (result) 


C*2  (»+3)/2.2-B 
rmis {initial  sequence) 


(3.5) 


where:  rms (error)  is  the  rms  value  of  the  computation 

error 

rms (result)  is  the  rms  value  of  the  resulting 
spectrum 

C is  a constant  that  depends  on  the 
machine  arithmetic  and  varies 
between  0.4  and  0.9 
M is  the  number  of  stages  in  the  FFT 
(M  = Log2  N where  N is  the  number 
of  FFT  data  words,. so  N=  2M) 

B is  the  maximum  number  of  bits  in 
the  representation  of  a number  in 
the  machine  (word  length) 

rms (initial  sequence)  is  the  rms  value  of  the  time  series, 

scaled  so  that  the  magnitude  of  the 
maximum  time  series  value  is  equal 
to  unity.  (Thus  it  is  proportional 
to  the  inverse  of  a signal  "crest 
factor",) 

The  rms  error  in  dB  of  the  fixed  point  FFT  may  be  deter- 
mined from  the  signal  plus  noise  to  noise  ratio  as, 

?(M+3)/2  _-B 

FFT„ERROR(dB)  <20  Log  (1  + - ■ n — rn = . ( 

F rms  (initial  sequence) 
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For  a transient  signal  f (t) , rms (initial  sequence)  in 


this  context  is; 


rms (initial  sequence)  = 


f*(t)  dt 


fmax(t)>  T 


or  in  terms  of  acoustic  parameters 

1 (Total  "intensity")  ^ 

T Peak  ‘‘signal" 

where  T is  the  duration  of  the  observation  time  interval. 

For  shallow  SUS  shots  rms (initial  sequence)  is  of  the  order 
of  0.003/T,  unfiltered  (0.01/T,  filtered  at  1 kHz).  Thus  the 
fixed  point  FFT  error  will  be; 

c#2((M+3)/2)-B 


FFT^ERROR (dB ) <20  Log  |1  + 


3 • 10~3/T 


(3.7) 


,12, 


For  a 16  bit  word  length,  a 4096  (2  ) word  FFT  of  1 second 

duration,  and  a 1 kHz  filter,  this  yields  a maximum  error  of, 

n ((12+3)/2)~16 
FFTpERROR  < 20  Log  jl  + u’--  - — OI7f 

= 20  Log  1.04 
= 1 dB 


A lower  bound  estimate  for  the  error  of  a spectrum,  from  Ref 
erence  10  gives, 

FFTpERROR  (dB)  > 20  Log  (1  + (M-2.5)  **  ( . 3)  2'B) 

which,  for  the  above  example,  gives  a minimum  rms  error  of 
the  order  of  0.005  dB. 

For  the  experimental  data  given  in  the  same  paper  from 
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Reference  10,  the  error  was  between  one  tenth  of  the  upper 
bound  and  the  upper  bound. 

The  error  analysis  expression  predicts  increased  error 
for  increased  signal  "crest  factor",  decreased  machine  word 
length,  and  increased  number  of  PFT  words. 

For  a floating  point  FFT,  upper  bound  of  the  error  is 
typically  less,  and  is  given  byj 

FFTpiERROR(dB)  = 20  Log  (1  + 0. 5‘2"ra»MJs) 

where  m is  the  number  of  bits  used  in  expressing  the  mantissa 
of  the  floating  point  word,  and  M as  before  is  the  number  of 
stages.  For  comparison  with  the  previous  example,  a machine 
would  use  two  integer  words  to  represent  a single  precision 
floating  point  word,  so  that  perhaps  20  bits  might  be  used 
in  the  mantissa,  giving  an  error  of  about  3»10~^  dB,  which 
is  completely  negligible  for  this  analysis.  The  increased 
precision  of  the  floating  point  FFT  is  achieved  at  the  expense 
of  substantially  increased  machine  memory  (2x)  and  a much 
longer  computation  time. 

Figures  3-8a  and  b give  plots  comparing  the  FFT  and 
exact  spectra  for  a set  of  conditions  that  correspond  to 
the  above  examples,  performed  on  a Varian  620L  minicomputer 
with  16  bit  arithmetic,  using  a fixed  point  FFT  routine. 

The  signal  analyzed  corresponds  to  a shallow  shot,  direct 
signal  only  (1.8  lb,  60  ft  depth,  300  ft  range)  1024  Hz  cut- 
off frequency,  samplinq  rate  4096  Hz  (fg  = 2 x FN) , 1 second 

36 


of  data,  16  bit  representative) . The  agreement  between  the 
FFT  and  true  spectrum  is  so  good  that  these  two  curves  super- 
impose on  each  other.  The  dashed  line  shows  their  differences 
on  an  expanded  scale  while  the  triangles  show  the  differences 
in  the  1/3  octave  band  spectral  levels.  Computations  on  the 
difference  data  indicate  an  rms  error  from  1 to  500  Hz  of 
about  0.3  dB  (excluding  the  0 Hz  term) , which  compares  favorably 
with  the  above  fixed  point  FFT  upper  bound  error  estimate  of 
1 dB  . The  mean  bias  of  about  0.1  dB  in  both  the  1 Hz  spec- 
trum and  the  difference  in  1/3  octave  band  values  is  believed 
to  be  the  result  of  round-off  error  in  the  FFT. 


3.2.3  Aliasing  Error.  Figures  3- 9a  and  b show  the  effect  of 
increasing  the  antialiasing  filter  cutoff  frequency  to  2048  Hz, 
which  is  also  the  Nyquist  frequency  fo^  this  sample  rate.  The 
positive  bias  introduced  is  clearly  obvious  at  high  frequency 
(above  250  Hz) . Using  the  theory  of  Appendix  B the  predicted 
error  (ratio  of  signal  plus  noire  to  noire)  is: 


El 1 * (f) 

ALIAS  ERROR(dB)  = 10  Log 


fc  2“  C 
= 10  Log  1 + — ar— 

rN  N 


f < f « f , f < 
r c C — N 


(3.8) 


For  f = 1024  Hz  (i.e.  fg  > 2 fN)  , Equation  3.8  predicts 
an  alias  error  of  less  than  6»10  ^ dB,  which  cannot  be  detected 
in  Figure  3-8  in  the  presence  of  the  FFT  round-off  error. 
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fjj  = 2 kHz,  Equation  3.8  predicts  an 


For  f = 500  Hz,  f = 

c 

aliasing  error  of  about  0.2  dB,  which  is  comparable  to  a 0.1  dB 
increase  in  error  observed  in  the  500  Hz  1/3  octave  band 
(compare  Figures  3-8b  and  3-9b) . The  increase  in  negative 
going  errors  at  frequencies  below  250  Hz  in  Figure  3-9b  is 
attributed  to  interaction  with  the  round-off  error  term. 

3.2.4  Finite  word  length  quantization  error.  Traditionally, 
when  dealing  with  noise  signals,  quantization  error  has  been 
treated  as  an  additional  noise  process  because  the  quantization 
process  can  be  represented  as  the  addition  of  a small  oscil- 
lating signal  (magnitude  equal  to  half  the  quantization  inter- 
val) that  has  a uniformly-distributed  amplitude  independent 
of  the  quantized  time  series.  For  a SUS  signal,  error  is  cor- 
related with  the  signal.  Here  we  find  that,  for  times  where 
the  signal  is  rapidly  varying  (i.e.  has  high  frequency  content) , 
it  also  has  a relatively  high  amplitude,  and  conversely,  low 
amplitude  is  correlated  with  low  frequencies,  so  the  signal 
cannot  be  uniformly  well-described  over  the  frequency  range  of 
interest  by  a few  discrete  steps.  For  the  shot  studied  (1.8  lb, 
60  ft  depth,  300  ft  range)  the  peak  pressure  (when  filtered 
to  1 kHz)  is  about  16  psi,  while  the  minimum  negative-going 
pressure  is  about  0.4  psi.  Thus  a 4- bit  representation  would 
completely  eliminate  (set  to  zero)  all  time  series  values 
except  those  in  the  immediate  vicinity  of  the  shock  wave  and 
bubble  pulse  peaks.  Since  the  long  periods  of  low  pressure 
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with  very  slow  pressure  gradients  carry  the  low  frequency 
information,  it  is  this  region  that  is  affected  most  by  low 
quantization  resolution. 

Figures  3-10a  through  3-13b  show  the  effects  of  progressive 
reduction  in  the  number  of  bits  (from  14  to  8)  used  to  represent 
the  signal  (see  Figure  3-8  for  the  16  bit  case).  (Note:  more 
than  50%  overrange  for  the  input  signal  was  allowed  (i.e.  peak, 
signal  less  than  66%  of  full  scale).}  The  error  shown  in  these 
figures  is  a combination  of  the  FFT  round-off  error  and  the  quanti- 
zation error.  Since  the  errors  are  of  a related  type,  there 
is  considerable  interaction  between  the  two.  Table  3rl  shows 
a comparison  of  single  number  measures  of  the  errors.  The 
trend  of  the  data  is  for  a slight  reduction  in  error  going  from 
the  16-bit  (Figures  3-8a  and  b)  to  14-bit  representation,  fol- 
lowed by  a progressive  increase  in  error  for  further  reduction 
in  number  of  bits.  The  FFT  input  data  were  not  corrected  to 
zero  mean  value  in  this  analysis,  so  that  the  substantial 
growth  in  the  DC  (0  Hz)  term  in  the  FFT'  is  also  indicative  of 
the  fact  that  quantization  is  not  an  unbiased  operator  for  this 
signal. 

The  net  decrease  in  error  of  the  FFT  for  a change  from 
16-  to  14-bit  signal  representation  is  interpreted  as  a reduc- 
tion in  FFT  round-off  error  while  still  preserving  adequate  signal 
representation.  Unfortunately,  the  FFT  error  analyses  in  the 


(Text  continues  on  Page  52) 


Figure  3-10a.  Effect  of  quantization.  Standard  shot,  all  parameters 
same  as  Figure  3-8  except  quantization  14  bits. 
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Figure  3-8  except  quantization  10  bits 
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literature  do  not  incorporate  this,  but  rather  deal  only  in 
fully-scaled  signals,  so  we  are  unable  to  predict  this  effect 
quantitatively.  It  can  be  understood  intuitively  from  the 
fact  that,  with  fewer  bits  in  the  incoming  signal,  there  is 
less  need  for  rounding  the  results  of  additions  in  the  early 
stages  of  the  FFT,  thus  reducing  the  errors  that  are  propagated 
to  later  stages. 

Table  3rl  gives  the  rms  error  for  1 Hz  band  spectra  and 
mean  error  in  1/3  octave  bands  for  the  calculations  shown  in 
the  previous  figures.  For  a 16-bit,  fixed  point  FFT,  a 14-bit 
signal  representation  gives  a minimum  in  combined  round-off  and 
quantization  error  and  a 12-bit  representation  is  almost  as 
precise  as  a 16-bit  representation.  Neither  of  these  conclusions 
is  expected  to  apply  to  a floating  point  FFT  since  that  would 
have  a much  lower  FFT  round-off  error.  In  the  floating  point 
case,  increasing  the  number  of  bits  in  the  signal  representation 
is  expected  to  provide  a more  precise  FFT;  the  rms  errors  in 
Table  3-1  would  certainly  be  upper  bounds.  Further,  the  errors 
for  1/3  octave-  data  are  expected  to  represent  accurate  estimates 
of  mean  error  for  a floating  point  FFT,  since  the  indicated 
biases  are  in  large  part  due  to  data  quantization,  and  so  are 
not  strongly  dependent  on  the  type  of  FFT. 

This  source  of  error  can  apparently  be  made  sufficiently 
small,  using  readily  available  hardware,  to  be  negligible  for 
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the  purposes  of  this  study  (e.g.,  minimum  12  bit  "sample  and 
hold"  A-D  with  minimum  4 kHz  sample  rate) . 

We  further  note  an  interaction  error  that  is  present  on 
some  A-D  converters/  between  sampling  and  quantization  errors 
which  perform  the  A-D  by  comparing  against  a "running"  signal. 

The  error  manifests  itself  in  the  form  of  time  base  "jitter", 
due  to  the  fact  that  this  type  of  converter  compares  against 
the  continuously-changing  analog  signal.  Thus  the  particular 
time  instant  at  which  the  sample  value  is  taken  is  unknown, 
except  that  it  is  bracketed  into  a time  window  of  some  speci- 
fied duration.  This  problem  is  completely  avoidable,  through 
hardware  selection  that  uses  “sample  and  hold"  circuitry,  and 
was  therefore  not  considered  in  the  analysis. 

3.2.5  Registration  effects  and  other  signal  processing  questions. 
Registration  is  a term  used  in  the  context  of  signal  proces- 
sing to  refer  to  the  temporal  location  of  the  signal  within 
the  time  window  with  respect  to  the  time  of  some  other  event, 
such  as  the  analysis  start  time  or  the  time  of  occurrence 
of  another  signal.  Since  the  analysis  of  the  DFT  processes 
as  applied  to  transient  events  (Section  3.1)  does  not  indicate 
the  presence  of  any  dependence  on  timing  of  the  signal,  other 
than  that  the  signal  is  completely  contained  within  the  time 
window  for  the  analysis,  we  expect  no  dependence  on  timing. 
Further,  since  the  spectrum  results  from  an  error-free  FFT 
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(no  alias,  FFT  round-off  or  quantization  errors)  that  is  exactly 
correct  and  interpolatable  exactly,  any  time  domain  information 
with  respect  to  two  or  more  signals  present  in  the  analysis 
window  can  be  recovered  exactly. 

When  error  sources  are  present  within  the  FFT,  registra- 
tion may  be  used  to  some  advantage  to  modify  these  errors. 
Figures  3-14  through  3-16  illustrate  these  effects.  In  Figures 
3-14  through  3-15,  the  signal  is  sampled  with  a h-  and  ^-sample- 
period  time  shift,  while  in  3-16a  the  signal  is  shifted  1,000 
samples  in  the  time  domain  by  pre-filling  the  data  sample  with 
zeros.  In  Figures  3-14  through  3-15  we  see  subtle  changes 
in  the  difference  spectrum,  but  the  mean  and  rms  errors  are 
not  significantly  altered  from  the  corresponding  values  for 
the  curve  in  Figures  3-8.  In  Figures  3-16  there  is  in  fact 
a statistically-signif icant  reduction  in  mean  and  rms  error, 
which  is  attributed  to  a reduction  in  s'FT  round-off  error  in 
the  early  stages  of  the  FFT  calculation  of  low  frequency  com- 
ponents (i.e.,  less  than  1 kHz). 

Two  other  aspects  of  digital  signal  processing  are  appro- 
priate to  mention  at  this  time.  The  first  concerns  the  addition 
of  zeros  to  round  out  the  length  of  the  time  series  to  some 
convenient  value,  while  the  second  concerns  the  use  of  windowing 
functions.  Equation  3.1  clearly  indicates  that  the  addition 
or  deletion  of  terms  for  which  p(qAt)  has  zero  value  has  no 
effect  whatsoever  on  the  value  of  the  resulting  spectrum. 

It  is  precisely  this  lack  of  effect,  coupled  with  the  finite 
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time  duration  of  the  transient,  that  permits  the  exact  reali- 
zation of  the  summation  with  infinite  end  points.  Further, 
when  the  time  series  is  fully  contained  within  the  analysis 
time  window,  there  are  no  signal  discontinuities  at  the  end 
points  that  are  not  already  part  of  the  signal,  and  therefore 
there  is  no  need  to  "smooth"  the  signal  by  applying  a window 
such  as  a Hamming  or  hanning  smoothing  function.  Since  the 
signal  is  by  definition  non-stationary , the  application  of 
such  a weighting  must  give  a result  that  depends  on  where  the 
signal  is  in  the  window,  with  no  hope  of  being  able  to  appro- 
priately correct  the  result  for  the  weighting. 

3.3  Surface  Reflections 

Having  used  the  analytic  model  to  examine  signal  proc- 
essing parameters,  it  is  appropriate  to  ask  what  other  information 
might  be  learned.  Two  effects,  clipping  of  the  surface  reflec- 
tion and  change  of  the  spectrum  shape  with  range,  have  been 
examined  and  are  discussed  in  this  and  the  subsequent  subsection. 

The  effect  of  surface  reflections  on  the  resulting  data 
is  of  critical  importance  for  shallow  (60  ft)  shots  because 
there  is  no  receiver  location  for  which  the  direct  signal 
can  be  isolated  in  the  time  domain  from  the  reflected  signal. 

A further  apparently  unavoidable  problem  of  shallow  shots  is 
the  fact  that  the  surface-reflected  shock  wave  frequently  is 
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truncated  because  of  cavitation  and  the  total  signal  may  be 
further  modified  by  the  non-planar  sea  surface.  In  the  model, 
because  of  lack  of  detailed  data,  this  is  treated  as  a truncation 
("clipping")  of  the  shock  wave  at  some  constant  value  with  no 
compensation  to  the  total  impulse. 

When  the  surface-reflected  signal  is  an  exact  replica  of  the 
direct  signal  with  the  exception  of  a fr  phase  change  and  a time 
delay  t,  the  resulting  signal  may  be  represented  as  follows: 


ptotal^  ~ pdirect^  + psur  ref^ 

* Pdireot(t)  ' Pdireot(t"T) 

= Paired  <1  - e+i“T> 
the  energy  spectrum  takes  the  form; 

Etotal<“>  * Edirect<“>  ‘2  - <e+i“T  + e'iWT»> 

* Edirect(“>  - cos  ur)> 

' 4 Edir.ct<“>  si"2 


Thus  the  sum  of  direct  and  linearly  surface-reflected  signals 
gives  an  energy  spectrum  which  is  modulated  by  a term  of  the  form 
4 sin  (^-).  This  term  has  alternate  maxima  and  minima  at  fre- 
quencies corresponding  to  ujt  = nv , n=0,l,2.  . . At  the  maxima 
(odd  values  of  n) , the  total  energy  spectrum  is  6 dB  above  the 
direct  spectrum  alone,  while  at  the  minima  (even  or  zero  n) , the 
interference  is  complete.  Thus  the  linear  reflection  case  is 
well  defined  and  accurately  calculable  from  frequency  domain  infor- 
mation, except  in  the  frequency  intervals  in  the  vicinity  of 
minima  of  the  sin  (sp-)  term.  The  principal  problems  for  data 
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analysis  are  thus  expected  to  be  those  associated  with  non-linear 
reflection  processes  (i.e.  signal  clipping). 

To  illustrate  the  effect  of  this  clipping,  the  combined 
exact  signal  spectrum  for  the  direct  plus  clipped,  reflected 
signal  is  compared  with  the  corresponding  unclipped  exact 
spectrum.  Figures  3-17  through  3-21  show  these  comparisons 
for  the  case  of  the  standard  shot  (1.8  lb,  60  ft  depth,  300  ft 
range,  shock  peak  pressure  42  psi)  combined  with  a surface 
reflection  that  has  no  attenuation  relative  to  the  direct 
signal  except  for  a truncation  of  the  shock  wave  to  15  psi 
amplitude.  The  delay  t between  the  direct  and  reflected  signal 
is  the  variable  in  this  group,  varying  from  24  ms  to  1.5  ms, 
corresponding  to  far-field  radiation  angles  of  90°  to  7° 
respectively.  As  can  be  clearly  seen,  the  sin  (j^)  modula- 
tion dominates  the  general  shape. 

Of  considerable  interest  is  the  fact  that  the  clipped  signal 

2 

spectrum  is  significantly  altered  only  at  the  minima  in  the  sin 

modulation  with  relatively  little  change  in  wide  band  energy  at 

other  frequencies.  The  fact  that  the  spectral  changes  due  to 

clipping  occur  primarily  only  at  the  spectral  minima  is  shown 

by  the  relatively  small  change  (less  than  1 dB)  in  1/3  octave 

2 

band  levels  except  at  those  bands  containing  the  sin  modula- 
tion minima.  This  consistent  behavior  strongly  suggests  that 
the  detection  of  a "clipped"  surface  reflection  may  be  difficult 
using  only  spectral  information  in  a system  with  a noise  floor, 
but  also  that  the  effects  of  clipping  on  the  low  frequency  1/3 
octave  band  spectrum  may  not  be  as  large  as  previously  suspected. 


One  significant  departure  from  the  above  noted  is  the  20  Hz 
1/3  octave  band  for  short  delays,  where  the  error  exceeds 
3 dB  for  1.5  ms  delay.  No  conceptual  explanation  for  this 
unique  behavior  has  yet  been  identified. 

{Text  continued  on  Page  76) 
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Figure  3-i9a.  Same  as  3-17  except  delay  is 


3.4  Effect  of  Range 

As  noted  in  the  introduction  to  this  section,  one  of  the 
significant  factors  which  bears  on  the  usefulness  of  explo- 
sive signals  for  use  as  sources  in  propagation  experiments 
is  the  fact  that  the  shock  wave  portion  of  the  time  series 
signal  attenuates  more  rapidly  than  a small  amplitude  acoustic 
signal,  at  least  for  short  ranges.  Previous  authors  have 
experimentally  determined  the  rate  of  change  of  amplitude  and 
decay  constant  of  the  exponentially-decaying  shock  wave  vs 
range.  Since  these  experimentally-determined  terms  were  used 
in  the  analytic  model  of  the  signal  used  in  this  study,  the 
model  can  then  be  used  to  examine  the  way  in  which  this  model 
of  the  frequency  spectrum  of  the  signal  changes  with  range. 

In  particular,  we  seek  an  answer  to  the  question,  "If  the  shock 
wave  propagation  follows  the  empirically-determined  attenuation 
law,  what  is  the  magnitude  of  the  changes  in  the  spectrum  shape 
to  be  expected  with  range?"  The  answer  to  this  question  is 


I 


then  used  to  phrase  an  answer  to  the  question,  "What  is  the 

% 

expected  variation  from  spherical  spreading  of  the  source  level 
spectrum  as  a function  of  range?" 

Figures  3-22  and  3-23  show  calculated  shock  spectra  at 

i 

ranges  of  1 kyd  and  10  kyd  compared  to  the  baseline  spectra  at 
100  yd  for  the  direct  arrival  only  of  a 60  ft  depth  1.8  lb  shot. 

The  computations  were  made  using  the  exact  SPECTRUM  program  with 
data  scaled  back  to  100  yds  assuming  spherical  spreading.  The 
differences  between  the  spectra  reflect  only  the  dependence  on 
the  programmed  difference  of  the  time  series  versus  range. 

v 

These  differences  from  spherical  spreading,  as  may  be  noted  from 

Appendix  A,  are  a decrease  in  the  shock  wave  peak  in  proportion  to 
- 13 

(range)  * and  an  increase  m the  exponential  decay  time  constant 

1 0 22 

in  proportion  to  (range)  ’ . The  nominal  computed  peak  amplitudes 

of  the  shock  wave  at  these  ranges  (assuming  spherical  spreading) 

are  3.05  psi  and  0-226  psi  respectively.  With  the  exception  of 

the  slow  conversion  of  several  minima  in  the  100  yd  spectra  into 

minor  peaks  in  the  1000  and  10,000  yd  spectra,  this  seemingly 

large  change  in  the  time  series  wave  form  had  a remarkably 

subtle  effect  on  the  spectra. 

The  present  position  of  many  workers  in  the  field  is  that  a 
transition  from  shock  propagation  (i.e.  spreading  loss  not  pro- 
portional  to  (range)  ) to  acoustic  propagation  (spreading  loss  pro- 
portional to  range”^)  would  be  expected  to  occur  in  the  vicinity 
of  1 psi  peak  pressure  in  the  shock  wave.  One  type  of  experiment  to 

t % 

evaluate  the  presence  of  this  occurrence  would  involve  a comparison 

% 
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of  i/3  octave  spectral  levels  which  have  been  corrected  back 
to  a common  range  via  spherical  spreading.  Figure  3-24  pro- 
vides a plot  of  the  results  to  be  expected  from  an  error- free 
performance  of  such  an  experiment,  performed  under  the  assump- 
tion that  shock  wave  propagation  continues  to  a range  beyond 
that  shown  in  the  figure.  We  have  plotted  1/3  octave  energy 
source  level  spectrum  difference  from  that  at  100  yds,  with 
band  center  frequency  as  a parameter  in  this  figure.  The 
result  of  a transition  to  acoustic  propagation  on  such  a figure 
will  show  up  as  a transition  of  all  curves  to  a horizontal  slope 
at  some  common  range.  In  the  area  of  interest  (i.e.  1000  yd 
range  to  10,000  yd  range,  20  Hz  - 500  Hz) , the  slope  of  all 
curves  is  less  than  1 dB  per  decade  of  range  change.  This 
already  small  slope,  when  coupled  with  the  probable  magnitude 
of  experimental  errors  associated  with  propagation  loss  correc- 
tion, source  yield  variation,  etc,  strongly  suggests  that  this 
form  of  data  reduction  will  not  yield  a powerful  test  of  trans- 
ition range. 
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4.  MINI-EXPERIMENTS 


4.1  Introduction 

The  success  of  or  need  for  the  source  level  experiment 
is  contingent  upon  the  amount  of  knowledge  that  is  at  hand  during 
the  decision  and  planning  phases.  Such  items  as  the  sensitivity 
of  results  to  A-D  sampling  rates,  low  pass  filter  frequencies, 

1/3  octave  levels,  reflection  deconvolution,  and  measurement 
range  have  been  discussed  above.  The  repeatability  of  the 
acoustic  level  from  a SUS,  the  interaction  of  the  SUS  reflections 
at  the  surface  and  the  water  medium,  and  the  effect  of  measure- 
ment range  should  be  examined  further.  Although  modeling  of 
the  last  two  effects  was  undertaken,  the  results  are  only  esti- 
mates, since  the  exact  physics  controlling  these  phenomena  are 
poorly  known.  It  seems  appropriate  that  these  three  items  be 
examined  by  means  of  limited  scope  experiments,  mini-experiments, 
to  better  quantify  the  achievable  accuracy  in  source  level 
measurements. 


4.2  SUS  Repeatability 

The  repeatability  of  the  acoustic  level  from  a standard 
SUS  has  been  a question  mark  in  most  scientists*  minds  for  a 
long  time.  Such  variables  as  amount  of  explosive  detonated, 
burst  strength  of  the  case,  and  potential  aluminum  vaporization 
are  among  the  more  often  heard  problems  concerning  SUS  repeat- 
ability. 
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The  amount  of  explosive  detonated  is  a function  of  several 
things.  The  first  question  arises  during  casting  the  charge. 

The  full  1.8  lbs  is  not  cast  at  one  time  but  is  the  result  of 
two  or  three  separate  pours,  the  last  one  or  two  being  made 
to  make  up  for  the  shrinkage  of  the  TNT  when  it  solidifies. 

Thus  one  has  the  majority  of  the  block  as  one  unit  with  several 
thin  layers  poured  on  top.  It  is  a well  known  fact  in  explo- 
sive experience  that  if  a charge  is  cracked,  the  charge  repre- 
sented by  pieces  that  cracked  away  from  the  main  charge  may 
or  may  not  detonate.  If  they  do  detonate  it  can  be  with  a 
slight  time  delay  or  with  a low  order,  slow  velocity  detonation. 
The  separate  pours  for  the  charge  mav  result  in  a situation 
similar  to  the  cracked  charge.  There  is  also  the  possibility 
of  the  TNT  charge  cracking  between  the  time  of  the  pour  and 
the  time  it  is  used.  The  amount  of  explosive  actually  detonated 
affects  two  phenomena  - the  overall  level,  and  the  bubble  pulse 
frequency.  An  overall  chanqe  in  level  is  proportional  to  the 
cube  root  of  the  weight.  Thus  a change  from  1.8  to  1.6  lbs 
would  result  in  a 0.3  dB  change,  which  raav  or  mav  not  be  con- 
sidered tolerable.  A change  in  the  bubble  pulse  frequency  can 
be  more  serious  because  of  the  octave  processing  that  is  pres- 
ently in  vogue.  In  the  lower  1/3  octaves,  i,e.,  10  to  100  Hz, 
the  1/3  octaves  are  so  narrow  as  to  hardly  span  a complete  amp- 
litude cycle  in  the  frequency  domain.  Therefore,  the  shifting 
of  a null  into  or  out  of  the  fixed  1/3  octave  band  can  result 
in  significant  errors.  This  situation  is  similar  to  the  bubble 
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pulse  variation  as  a function  of  depth,  where  a few  percent 
change  in  depth  can  result  in  1/3  octave  level  changes  of  the 
order  of  0.5  to  3 dB. 

The  effect  of  the  burst  strength  of  the  case  on  the  acoustic 
source  level  has  been  questioned  in  recent  years.  The  situation 
can  be  paraphrased  as  follows:  energy  that  is  expended  to  burst 
the  case  containing  the  explosives  cannot  contribute  to  the 
acoustic  source  level.  Further,  we  know  that  there  is  a wide 
range  of  burst  strengths  for  the  case  and  this  will  cause  the 
acoustic  output  to  vary.  The  size  of  this  effect  is  difficult 
to  estimate,  but  in  any  case  this  concern  by  some  members  of 
the  community  should  be  borne  in  mind  when  planning  a source 
level  variability  experiment. 

Another  concern  that  is  heard  from  time  to  time  is  the 
possibility  of  vaporization  of  the  aluminum  case.  One  method 
of  enhancing  the  low  frequency  components  of  an  explosive 

I 

spectrum  is  to  add  uniformly-distributed  fine  aluminum  powder 
to  the  explosive.  The  result  is  a greatly  increased  bubble  due 
to  the  aluminum  vaporization  and  resultant  increase  in  bubble 
pulse  energy,  which  results  in  a much  higher  spectral  level  at 
the  lower  frequencies  than  would  otherwise  be  obtained. 

Following  this  line  of  reasoning,  if  the  TNT  i > the  SUS 
can  actually  vaporize  small  quantities  of  aluminum  from  the 
case,  the  source  level  at  the  lower  frequencies  could  be  affected. 
Another  comment  is  usually  added.  This  vaporization  will  take 
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place  only  part  of  the  time  to  add  to  the  uncertainty.  This 
concern  can  be  appreciated  and  is  a potential  cause  of  source 
level  variation,  but  is  difficult  to  evaluate. 

The  sources  of  small  variations  discussed  above  are  dif- 
ficult to  separate  except  in  very  elaborate  experiments.  The 
reasons  and  causes  of  shot-by-shot  variations  are  fine  to  con- 
template but  their  understanding  is  not  necessary  to  define  the 
repeatability  of  the  acoustic  source  level  of  a SUS  as  used  for 
propagation  loss  purposes.  Therefore,  it  is  proposed  that  a 
number  of  SUS  from  different  lots  be  detonated  and  their  relative 
outputs  measured  to  quantify  the  source  level  variation  without 
regard  for  potential  cause  of  the  error. 

4.3  Surface  Reflection  Non-Linearity 

When  shallow-detonated  SUS  are  used  as  acoustic  sources, 
the  surface  reflection  problem  is  ever  present.  Experiments 
purport  to  have  observed  the  reflected  shock  wave  truncated 
in  the  pressure  range  of  approximately  15  psi  to  essentially 
no  truncation,  as  is  theoretically  predicted  tor  deaerated 
water.  Such  a variation  in  the  behavior  of  the  surface  reflection 
brings  up  several  problem  areas. 

1.  Can  deconvolution  techniques  adequately  handle  the 
range  of  variation  in  the  surface  reflection? 

2.  How  does  this  variation  reflect  in  the  use  of  source 
levels  for  PL  measurements? 
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3,  Is  there  any  correlation  between  angle  of  encounter 

with  the  surface  and  the  shock  wave  truncation  process? 

The  truncation  process  for  the  reflected  shock  wave  can 
be  described  as  a cavitation  process?  i.e.,  when  the  negative 
pressure  of  the  shock  wave  overcomes  the  tensions  the  water 
can  support,  it  is  pulled  apart  with  the  forming  of  small  bubbles, 
or  in  the  case  of  large  charges  a layer  cavity  that  will  later 
collapse.  In  any  case  the  energy  that  is  represented  in  that 
part  of  the  shock  wa^e  that  has  been  ’•cutoff"  has  been  dissipated 
in  the  cavitation  process.  The  tensional  stress  that  water 
will  maintain  is  a function  of  the  dissolved  gas  in  the  water. 

For  totally  degasified  water  the  tensional  stress  that  is 
supportable  far  exceeds  the  tensional  stress  required  of  the 
medium  for  the  reflected  shock  wave.  However,  with  increasing 
gas  concent  the  tensional  stress  supportable  by  the  water 
medium  rapidly  decreases  so  that  truncation  of  the  shock  wave 
can  be  expected  as  a normal  occurrence.  Therefore,  the  expec- 
tation is  that  there  will  be  varying  truncation  of  the  shock 
wave  due  to  the  gas  content  of  the  water. 

Deconvolution  techniques  were  tested  in  the  analytical 
portion  of  this  study  and  the  results  presented  there.  One 
of  the  objectives  of  the  experiment  is  to  determine  whether  or 
not  the  deconvolution  error  estimates  developed  above  are  real- 
istic for  a real  measurement.  Such  a set  of  data  can  be  easily 
acquired  by  measuring  the  same  shot  at  two  locations,  one  where 


deconvolution  is  required  and  the  second  where  it  is  not,  i.e., 
the  reflection  arrives  after  the  second  bubble  pulse. 

The  second  area  of  concern  is  the  effect  of  variability 
in  the  surface  reflection  for  the  ray  paths  that  are  reflected 
from  the  surface  near  the  detonation  point.  The  amount  of 
variation  will  be  measured  so  that  an  estimate  of  the  variation 
ot  the  two  close-in  reflected  ray  paths  can  be  made  for  a 4-ray 
path  bundle. 

4.4  Effect  of  Measurement  Range 

A preferred  experimental  investigation  of  the  transition 
to  acoustic  propagation  would  be  self-correcting  for  propa- 
gation loss  and  shot  yield.  If  that  experiment  were  also  self- 
correcting  for  small  calibration  differences  in  the  transducer 
system,  the  accuracy  would  be  improved  even  further.  One  such 
experiment  would  involve  the  use  of  two  wideband  transducers 
to  receive  signals  from  shots  dropped  in  an  annular  region 
above  one  of  the  transducers,  which  is  placed  deep.  The  second 
transducer  is  placed  shallow  and  at  a significant  horizontal 
range  from  the  source  drop  area.  This  latter  range,  plus  the 
radius  of  the  drop  area,  is  selected  to  span  the  horizontal  range 
region  of  interest.  The  deep  transducer  is  used  to  record 
signals  at  nominally  constant  source  range,  as  a reference. 

To  illustrate  the  data  reduction  scheme,  we  require  an  estimate 
of  the  signal,  which  we  take  as  follows  (from  the  analytic  model 
used  here;  see  Appendix  A). 
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Shock  wave  peak  pressure  = P^  = 2.08  x 10^  (W1//3/R) 1,13 

3 1/3 

First  bubble  pulse  peak  pressure  = = 3,3  x 10  (W  ' /R) 


Using  our  analytic  model  (SHOTTIME)  as  a representation 
of  a band-limited  signal,  it  has  been  determined  that  these 
signals  for  a given  shot  can  be  detected  with  a relative  in- 
accuracy of  less  than  5%  for  systems,  with  an  upper  cutoff  fre- 
quency in  the  range  of  10-20  kHz.  The  need  for  such  a wideband 
system  arises  from  the  desire  to  study  peak  pressure  independent 
of  the  decay  rate  of  the  exponential  terms,  which  are  integrated 
into  the  data  in  narrower  bandwidth  systems.  The  data  reduc- 
tion consists  of  determining  the  peak  pressure  of  the  shock  wave 
and  bubble  pulse  at  each  of  the  transducers,  and  forming  their 
ratio.  The  result  at  each  transducer  if  shock  wave  propagation 
holds  is: 


P1(R1)/P2(R1) 


2,08  x 104 
3.3  x 103 


— 

(W1/J/R1) 


1.13 


= 6.3 


We  note  that  this  ratio  will  be  independent  of  transducer 
calibration  constant,  as  long  as  the  system  response  is  linear. 
We  r.ow  form  the  ratio  of  this  quantity  for  the  same  shot,  at 
two  different  ranges,  R^  and  RQ,  i.e.j 


P1(R1)/P2(R1}  = / w1/3/R! 

P1<R0)/P2(V  y W1/3/R0 

for  shock  wave  propagation  at  R^,  Rj. 


0.13 
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The  result  of  the  data  reduction  is  an  experimental  value 
that,  for  the  case  of  shock  wave  propagation,  is  now  inde- 
pendent of  shot  yield,  and  depends  only  on  the  relative  range 
of  the  two  measurements.  Consider  now  a case  where  both  and 
Rq  are  such  that  acoustic  propagation  of  the  shock  wave  occurs. 
One  representation  of  this  for  the  shock  wave  is; 


P1(R1)  = 2.08  x 104  (W1/3/Rm)1,13  (Rm/R) , R>Rm. 

This  expression  equals  shock  wave  amplitude  at  range  Rm, 
but  follows  acoustic  propagation  for  R > R^.  If  both  R^  and 
Rq  above  are  greater  than  R^,  then  the  ratio  is; 


P1(R1)/P2(R1) 

p-(R0)/p2(R0) 


(wl/3/v 

{wl</3/Rm) 


0.13 

7T.13 


l; 


rx,  RqsR 


m 


If  Rq  ' < R^  we  have; 


P1(R1)/P2(R1) 


(wl/3/Rm} 


0.13 


W/P2(R0J 


T7T;—' or 

0; 


<W~'J/Rn) 


0.13 

1;  VVR1 


In  the  above,  if  RQ  is  the  nominally  constant  range  distance, 
the  ratio  of  interest  is  nominally  constant  at  a value  less  than 
unity,  which  permits  computation  of  the  range  R^  whore  the 
transition  to  acoustic  propagation  occurs.  A plot  of  this 
data  vs  log  (variable  range  R^)  that  spans  the  transition  range 
R is  expected  to  show  a slope  of  -0.13  below  R . and  approach 

Cl  m 


a constant  value  above  R , 

m 


In  summary,  the  above  discussion  suggests  that  an  experi- 
ment to  evaluate  range  dependence  in  the  frequency  domain  is 
fraught  with  difficulty  because  of  the  small  differences  that 
are  under  study.  A direct  experiment  in  the  time  domain  is 
suggested,  and  a data  reduction  scheme  is  outlined  that  is  self- 
correcting  for  yield  and  range  and  that  could  provide  the  desired 
information.  Further  study  to  define  equipment  requirements 
is  necessary,  however,  before  a final  experiment  plan  can  be 
defined. 

4.5  Specific  Experiments 

During  October  1977,  PAR  wet  tests  were  performed.  To  take 
advantage  of  that  cpportunity,  an  experiment  plan  was  prepared 
to  investigate  SUS  repeatability,  decorrelation  and  non-linear 
surface  reflection  effects  and,  to  a lesser  extent,  the  effect 
of  measurement  range.  This  experiment  plan  is  given  in  Appendix 
C. 
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5.  CONCLUSIONS  AND  RECOMMENDATIONS 


The  foregoing  analyses  and  discussions  support  a number 
of  conclusions  concerning  the  errors  to  be  expected  in  an 
experimental  determination  of  SUS  Source  level. 

• Transducer  calibration  errors  of  less  than 
0.5  dB  can  be  achieved. 

• Signal  processing  errors  using  an  FFT  analysis 
system  can  be  reduced  to  less  than  0.1  dB  bias  in  1/3  oc- 
tave bands,  and  less  than  0.5  dB  rms  error  in  1 Hz  bands  in  the 
frequency  range  12.5  Hz  to  500  Hz,  using  a system  consist- 
ing of: 

1)  Minimum  12  bit  or  more  Analog  to  Digital  con- 
verter with  negligible  time  base  jitter  (to  min- 
imize quantization  error) . 

2)  Antialias  filtering  at  a cutoff  frequency 

of  one-half  or  less  of  the  Nyquist  frequency  with 
at  least  a 24  dB/octave  roll-off  of  the  low-pass 
filter  skirt  (to  minimize  aliasing  error) . 

3)  Rectangular  windowing  of  the  time  series  with 
the  entire  event  included  in  the  window. 

4)  An  FFT  algorithm  that  uses  at  least  16  bit 
fixed  point  computation  representation  (to  minimize 
computation  roundoff  error) . 

• Physical  knowledge  of  the  laws  governing  a transition 
from  shock  to  acoustic  propagation  is  presently  limited, 
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but  even  shock  propagation  produces  less  than  a +0.3  dB 
range  dependence  for  the  1/3  octave  spectrum,  in  the 
range  of  1000  - 10,000  yds.  Taking  this  as  an  upper- 
bound  to  the  variability  due  to  range  dependence  provides 
one  estimate  for  this  source  of  variability  that  should 
be  verified  through  a range  dependence  mini-experiment. 

• Similarly,  physical  understanding  of  the  variability 
of  spectra  due  to  the  deconvolution  of  a signal  including 
a direct  signal  and  a clipped  surface-reflected  signal 
suggests  that  while  the  variations  in  1 Hz  bands  are 
large,  the  1/3  octave  band  spectra  vary  from  the  idealized 
(i.e.,  unclipped)  spectra  by  less  than  0.5  dB  on  the  average 
in  the  frequency  range  25  Hz  - 250  Hz  without  any  cor- 
rection for  clipping  (see  Figures  3-16b  to  3-2Qb) . This 
suggests  strongly  that  some  very  simple  analysis  algo- 
rithms based  on  propagation  delay  only  can  be  used  to 
“deconvolve"  shallow  shot  data  regardless  of  clipping. 

This  suggestion  must  also  be  tested  in  the  physical  world, 
and  a mini-experiment  is  described  in  this  report  to  eval- 
uate it. 

Using  the  above  error  estimates  permits  an  estimate  of 
overall  confidence  level  in  an  experimental  determination  of 
SUS  source  levels. 

We  assume  that  the  above  errors  are  random  and  non-inter- 
acting  and  sum  their  squared  values  to  arrive  at  a mean  square 
error  estimate  for  the  completed  experiment.  One  additional  fac- 
tor that  deserves  consideration  is  the  data  record-playback 
system. 
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We  estimate  that  with  current  state  of  the  art  and  individual 
gain  setting  calibrations  that  this  error  should  be  less  than 
0.5  dB.  The  root-sum-squared  of  the  above  five  errors  (range 
0.3  dB,  calibration  0.5  dB,  record-playback  0.5  dB,  processing 
0.1  dB,  deconvolution  with  clipping  0.5  dB)  i3  about  0.9  dB. 

This  indicates  to  us  the  potential  feasibility  of  an  experi- 
mental determination  of  SUS  source  level  to  within  an  absolute 
error  of  +1  dB,  subject  only  to  the  positive  verification  of 
the  error  limits  for  surface  clipping  and  range  dependence  in 
mini-experiments . 

The  above  error  estimate  applies  to  the  certainty  with 
which  the  source  level  from  a carefully  observed  SUS  shot  can 
be  known.  One  further  factor  in  the  uncertainty  in  using  this 
data  to  predict  the  source  level  of  another  shot  is  the  yield- 
depth  variability  from  shot  to  shot.  A mini-experiment  to  shed 
further  light  on  these  statistics  for  SUS  is  also  described. 

Since  there  is  no  data  or  a priori  evidence  currently  available 
on  this  portion  of  the  uncertainty,  experimental  study  is  our 
only  approach  to  this  factor.  Since  it  could  have  some  influence 
on  future  usability  of  a SUS  calibration,  it  should  be  undertaken 
prior  to  final  commitment  to  the  calibration  experiment. 

Thus  our  recommendations  from  this  study  concerning  SUS 
calibration  include; 

• Perform  a mini-experiment  to  identify  bounds  for 

range  dependence  of  SUS  spectra. 

• Perform  a mini-experiment  to  identify  the  bounds  for 
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error  due  to  clipping  and  deconvolution  of  shallow  shot 
spectra. 

• Perform  a mini-experiment  to  evaluate  the  statistics  of 
shot-to-shot  variability  of  SUS. 

• Evaluate  the  significance  of  the  data  from  the  above 
experiments  prior  to  committing  to  the  performance  of  an 
experiment  to  improve  SUS  calibration  accuracy. 
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APPENDIX  A 


AN  ANALYTIC  FUNCTION  MODEL  OF  AN  EXPLOSION  SIGNAL 
AND  PROGRAMS  FOR  COMPUTATION  OF  ITS 
TIME  SERIES  AND  SPECTRUM 


APPENDIX  A 


Ah  Analytic  Function  Model  of  an  Explosion  Signal  and 
Programs  for  Computation  of  its  Time  Series  and  Spectrum. 

The  explosion  signal  used  in  this  model  is  intended  to 

represent  a shock  wave  and  four  successive  bubble  pulses.  The 

model  uses  all  currently  available  scaling  information  for 

describing  the  time  series  signal. 

The  characteristics  included  in  the  signal  representation 
2 

are  the  following: 

• peak  pressure  of  shock  wave 

0 peak  pressure  of  four  bubble  pulses 

0 maximum  first  negative  phase  pressure 

0 time  of  zero  crossings  during  shock  wave 

and  first  bubble  pulse. 

0 exponential  decay  constant  of  shock  wave 
0 time  to  peaks  of  first  four  bubble  pulses 
The  decision  to  represent  the  time  series  as  piecewise 
analytic  using  simple  functional  forms  was  based  on  the  desire 
to  make  an  exact  Fourier  transform  of  the  time  series  so  that 
the  spectrum  could  be  known  analytically,  and  further,  that 
Fourier  transform  techniques  could  be  used  to  provide  a simu- 
lation of  anti-alias  filtering . A similar  approach  for  signal 
modeling  was  taken  in  Reference  11  independent  from  this  effort. 

2 ‘ 

J.  B.  Gaspin  and  V.  K.  Shuler,  op.cit. 
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A. 1 The  Unfiltered  Time  Series 

Figure  A-l  shows  the  components  of  the  signal  model.  The 
time  axis  was  divided  (at  the  cusps  in  the  time  series  curve) 
into  epochs,  and  four  different  functions  defined  in  each  epoch 
are  summed  to  define  the  time  history.  The  four  functions 
consist  of  an  exponential  decay  at  the  beginning  of  the  epoch, 
an  exponential  rise  at  the  end,  a positive  constant  and  a 
negative-going  half  sine  wave. 

In  addition  to  the  peak  pressure  and  time  values  identi- 
fied earlier,  several  additional  constraints  and  assumptions 
were  imposed  in  order  to  define  the  values  of  the  constants, 
including; 

1.  Pressure  is  continuous  across  the  end  of  each  epoch 
(a  physical  constraint) . 

2.  Zero  mean  value  of  pressure  was  assumed  within  each 
epoch  (except  that  epoch  five  is  merged  with  epoch 
four  to  apply  this) . 

3.  The  exponential  decay  constant  was  assumed  to  have 
the  same  value  on  both  sides  of  each  bubble  pulse. 

These  twelve  constraints  and  assumptions  plus  the  values 
of  the  thirteen  constants  supplied  as  data  v.’ere  still  insuf- 
ficient to  define  the  thirty  independent  constants  of  the  model, 
so  the  following  arbitrary  selections  of  values  were  made: 

PC(1)  « PC (2) 

PC  (3)  = PC  (4)  = 0 

6B  (2)  = 0B  (3)  = 6 B ( 4 ) = 0A<3)  0A(4)  = 0A(5) 


Components  in  the  time  domain: 

PA,  » PA(i)  expl (t (i) -t) /0A(i) ] 

PB,  » PB(i)  exp( (t-t  (i) /OB  (i) ] 

PC.  « PC(i) 

PD.  = PD(i)  ainl  (t-t  (i) ) nT(i)  ) 

i is  the  epoch  number 

t is  elapsed  time 

t(i  + l)  is  the  elapsed  time  to  the  end  of  epoch  i 

T(i)  is  the  duration  of  epoch  i. 

All  functions  are  defined  to  be  zeros  outside  their  epoch 
of  definition. 


Figure  A~1  The  Unfiltered  Time  Series 
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As  a result  of  these  considerations,  all  of  the  constants 
of  the  model  could  then  be  put  in  the  form  of  shot  parameter 
dependent  variables  through  the  similarity  expressions,  giving 
the  expressions  shown  in  Table  A-l. 

A. 2 The  Fourier  Spectrum  and  the  SPECTRUM  Program 

The  Fourier  spectrum  is  computed  from  the  standard  Fourier 

integral: 

S (w)  = / p(t)  e“-*wtdt  (A.  1 ) 

— Oo 

The  expressions  given  in  Section  A.l  were  used  for  p(t), 
with  tho  approximation  that  the  domains  of  the  exponential 
terms  were  extended  to  plus  or  minus  infinity  as  appropriate. 

This  approximation  is  believed  to  have  negligible  impact  on 
the  spectrum  so  lonq  as  the  exponential  time  constant  is  less 
than  about  one  fifteenth  of  the  bubble  pulse  interval.  This 
approximation  limits  the  applicability  of  the  program  to  shot 
parameters  such  chat; 

2 < 1700  ¥l/4 

For  a 60  ft.  SUS  shot  this  implies  a range  limit  _>f  about 
30,000  ft. 

The  seventeen  terms  of  p(t)  evaluated  analytically 

using  Equation  A.l.  l'he  results  for  the  generic  types  of  terms 
are  as  follows: 

Shock  wave: 

S_  (u>)  = PA ( 1 ) expUwtU)) 

s nmnT^ 


A~4 


U 


T 


Table  A-l 

Shot  Model  Parameters 

Z = D+33  where  D is  shot  depth  in  feet 
W = shot  yield,  lbs  of  TNT 
R » shot  range,  feet 
X =Wi/3/Z5/6 
Y » W1//3/R 
PI  = 3300  Y 


PA  ( 1 ) 

a 

2.08  x 104  Y1*13  - PC(1) 

PB  (1) 

a 

PI  - PC  €1) 

PC(1) 

a 

.403  PD ( 1) 

PD  (1) 

a 

8.0  Z2/3  Y 

0A{1) 

= 

5.8  x 10~5  Wl/3  Y"0*22 

0B(1) 

PB(l)”1  It  l - .4O3)PDU)T(l)-PA<li0A(l)] 

T (1) 

23 

4.34X,  t{2)  * 4.34X  ♦ t(l) 

PA  (2) 

S 

PB  (1) 

PB{2) 

« 

, 22P1  - PCU) 

PC  (2) 

S3 

PC(1) 

PD  (2) 

62 

2.42  (PCI  + PA(2)  exp(-. 135T(2)/0A{2) ) ) 

0A(2) 

SB 

0 B ( 1 ) 

OB  (2) 

as 

((_L_PD{2)  - PC  (2)  )T  (2)  - PA  ( 2 ) 0 A ( 2 ) ) /PB  ( 2 ) 

TT 

T(2) 

a 

3.06X,  t{3)  = 7.40X  ♦ t(l) 

PA  (3) 

83 

.22  PI 

PB(3) 

ta 

.10  PI 

A-5 
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Table  A-l  Con't 


PC (3)  » 0.0 

PD (3)  » .16  it  PI  6A(3)/T(3) 

0A{3)  = 0B(2) 

0B(3)  » 0B(2) 

T (3)  = 2. 48X,  t (4)  = 9.88X  + t(l) 

PA  (4)  = PB  (3) 

PB(4)  = Q.03  PI 

PC  (4)  « 0.0 

PD (4)  = 0.08tt  PI  AA(3)/T(4) 

0A (4)  = 0B(2) 

0B(4)  * 0B(2) 

T(4)  « 2.31X  t (5)  =*  12.19X  + t (1) 

PA (5)  = PB  (4) 

0A  (5)  = OB  (2 ) 
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Symmetric  bubble  pulse  of  amplitude  B and  decay  constant  6: 
S^w)  = 2B  exp(iu)t) 

e u/e*  + <w)2) 

Rectangular  wave  of  amplitude  C,  start  time  t,  duration  T: 
Sr(«)  * C exp(iait)  (exp(iu)T) -1) 


Half-sine  pulse  of  amplitude  D,  period  start  time  t, 
period  T: 

Sh(«)  = <-$-)  exP(i^)  (-exP 

( (-m~)  ~ w J 


A listing  of  an  interactive  Fortran  program  used  to  compute 
a 1000  point  spectrum  on  a Varian  620L  minicomputer  is  shown 
in  Table  A-2,  while  sample  output  listing  is  given  in  Table  A-3. 
The  program  will  compute  the  spectrum  of  a direct  signal,  or 
direct  plus  surface-reflected  signal  or  just  a surface -re fleeted 
signal,  with  an  option  to  truncate  the  magnitude  of  the  surface  - 
reflected  shock  wave. 

The  input  parameters  for  the  program  are: 

1.  Shot  yield,  depth  and  range  of  the  direct  signal 
(English  Units). 

2.  Reflected  signal  range,  time  delay  between  direct 
and  reflected  signal,  and  the  pressure  for  trunca- 
tion of  the  reflected  shock  wave. 

3.  Spectrum  parameters  (start  frequency,  number  of 
frequency  values  and  interval  between  frequencies. 
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The  output  listing  includes j 

1 . shot  data 

2.  a list  of  the  model  parameter  values 

3.  a listing  of  the  complex  Fourier  pressure 
spectrum,  the  energy  spectrum  and  the 
energy  level  spectrum 
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Table  A- 2 Listing  of  SHOT  SPECTRUM 


1 C NAME  SHOTSPECTRUK2A  **********RC= 32000  ************** 

2 C COMPUTES  FOURIER  SPECTRUM  OF  DIRECT  PLUS  REFLECTED 

3 C AND  OUTPUTS  RESULTS  TO  MAC  TAPE 

4 C 

5 C 

6 C IF  DIRECT  RANGE  ( RDX  =0.  DO  REFLECTED  ONLY 

? C IF  PCUP<=0,  DO  NOT  COMPUTE  REFLECTED 

8 C 

9 C STARTINC  FREQ  MUST  ALWAYS  BE  >=  0.0 

10  C START  TIHE=0  BY  DEFINITION 

H C 

12  c 

13  COMMON  T<  5) , TT(  4) , PA(  3) , PBl 4) , PC(  4) , PD(4) , AA(5) , AB(4) , XIMFPt 1) 

14  DIMENSION  TH( 5) ,REFP( 1000) 

15  COMPLEX  CFP, El. U2.U3.Z 

16  PI=3. 141592654 

17  EMAX*  0 . 

18  EI=(0. , l . ) 

19  C CEN*2(DYNE/SQ  CM/PSl)**2/R0C 

20  CEN=6 1720.368 

21  C INPUT  SHOT  DATA 

22  13  WRITE  (3.31) 

23  READ  (2,33)  W.D.RD 

24  IF< W.EO.0. ) CO  TO  1010 

25  WRITE  (3.32) 

26  READ  (2.33)  RR.TDL.CL1PP 

27  31  FORMAT( 43H  GIVE  YIELD(  LBS) ,DEPTH(  FT) .DIRECT  RANCE( FT) ./, 

28  *20H  SEPARATED  BY  COMMAS) 

29  32  FORMAT! 4511  GIVE  REF.  RANGE!  FT) .TIME  DELAY! SEC) .CLIP  PSI,/. 

30  '20H  SEPARATED  BY  COMMAS) 

31  33  FORMAT  (3F13.6) 

32  C INPUT  FREQ.  VALUES 

33  WRITE  (3.3) 

34  READ  ( 2 . 2) FO . FREQ, DFREQ 

35  2 FORMAT  (3V13.4) 

36  3 FORMAT! 56 H CIVE  FREQUENCY  PARAMETERS: START INC  FREQUENCY. NUMBER  OF 

37  *44HFRE0UENCIES  AND  INTERVAL  BETWEEN  FREQUENCIES/ l 5H  REAL. SEPARATED 

38  *32H  BY  COMMAS:  (MAX  512  FREQ  VALUES)) 

39  C PREPARE  OUTPUT  DEVICE 

40  WRITE  (3.4) 

41  4 FORMAT!  2411  MOUNT  DATA  TAPE  ON  MT0D 

42  PAUSE  10 

43  C PROCESS  INTERNAL  PARAMETERS  NEEDED  FOR  CALCULATION  MODEL 

44  C USING  TSVAL  SUBROUTINE 

43  T#30. 0 

46  DO  31  K' I. 1000 

47  REFPt K>  =0.0 

40  51  XINFP(K)  = 0.0 

49  C CHECK  TO  SEE  IF  REFLECTED  ONLY 

50  IF(  RD.l.K.0.0)  CO  TO  971 

5 1 PASS5 I . © 

52  R3RD 

53  TCC'0.0 

36  FCLIP'O.O 

35  TC=0.O 

36  61  CALL  TSV-Ul  W.D.R.T0) 

37  WRITE!  11.41) 

58  41  FORMAT! U3K  DIRECT  SHOT  PARAMETERS) 
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Table  A- 2 Con't. 


OUTPUT  DATA  HEADING 
WRITE! 11,5) 

WRITE!  1 1 , 7 > F0 . FREQ, DFREO 
WRITE!  1 1,40)T0 
WRITE!  1 1 ,6) W, D. R 
WRITE!  11,8) 

DO  0^0  ja  i , 4 

WRITE!  11*9) TV!  I)  ,PA!  I)  ,AA!  D.PB!  I),AB<  I)  ,PC(  I),PD(  I) 
WRITE  (1|, 10) PA! 5) , AA! 5) 

WRITE!  11.  12) 

DO  65  1=1.5 
TH(  I)s  I.xAA!  I) 


FIND  SPECTRUM  FROM  DATA 
COMPUTATION  START 
NST=  1 

0M=2.*Pl*F0 
OM0=OM 
NFREfl*  FREQ 
DELOM®  2 . *P l *DFR£Q 
CHECK  FOR  ZERO  FREQ 

953  1F10M.CT.0.0)  GO  TO  960 
COMPUTE  ZERO  FREQ 

RE3 PA!  1>*AA(  1)#F.XP(-TCC/AAC  l))+PC!  t) *AA( 3)*TCC*PCL1P 
RE=RE+PC!  1)*<T!3)-TC> 

RF.=  RE'PD!  1)*!  I . +COS! P I oTCCMT!  1 ) ) )/(  PI/TT!  1)  >+2.*PAl2)*AA(2> 
DO  954  K-2.4 
L-  K*  1 

954  RE«  RE>2 . *PB! K) * AA!  L) -2 . *PD<  K)  *TT< K) /P I 
REFPI NST) =REFP(  NST) ♦ PASS* RE 

OM=OM0*DE1.OM*FLOAT(NST)  , 

NST*NST*  1 ,i 

GO  TO  953 


FOR  POSITIVE  FREQUENCIES 
960  CONTINUE 

DO  9?0  I®NST, NFREQ 
Z3E1*0« 

CFP3 PA!  1 ) *CEXP! Z*TC-TCC*TR1  t ) ) /!  TH!  1 ) -Z) 

'♦PC!  1 ) *CKXP! Z*T<  3) ) /! TH!  3) -Z) 

'♦PC!  1>*<CEXP(Z*T!3)  >-CEXP!Z*TC> ) /Z 
'♦PCUP*CF.XPIZ*T!  l)  )*!  CEXP!'Z*TCC>-I.  )/Z 
XI“Pl/TT<  1) 

CFP=>CFP-PD!  l)*X(/!  X!**2.-0K**2. > *(  CEXP!  Z*T12) > ♦ 
'CF.XP!  7*TC ) *<  COS!  TCC*Xl  > -E l *S  I «( TCC*XI ) ) ) 

D!l  962  KM  . 4 
Ls  £♦  1 

962  CFP9CFP+PB! K) *U2( THi  L)  .OM.T!  L>) 

1>0  963  K®2.4 

l.°  K*  1 

963  CFP3C.FP-I‘»<  K)  *U3<  OM.T!  K>  .TT!  K) ) 

RKFP!  I ) “ REFIM  I > ♦REAL!  CFP)  'PASS 
XIMFp!  HsXIMFP!  1 ) ♦ A l MAC! CFP > *P  ASS 
EEST»RKrp«  I)»*2*XIMFP(  l)**2 

IF  ! ERST. CT. EMAX> EMAX* EEST 
Xl*t 

OMa  Xl*DEljOM»OMft 
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Table  A- 2 Con't. 


117 

11a 

119 

120 
121 
(22 

123 

124 
123 
526 
127 
126 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 

139 

140 

141 

142 

143 

144 
143 

146 

147 

148 

149 

150 

151 

152 
133 

154 

155 

136 

137 
(58 
139 
160 
161 
562 
H>3 
tot 

165 

166 

167 

168 

169 

170 

171 

172 

173 

174 


970  CONTINUE 

IF( PASS. EQ. -1 . ) GO  TO  1000 
C 
C 

C SECOND  PASS  NEEDED? 

C IF  NECESSARY.  PROCESS  CLIPPED  SIGNAL  INFORMATION 

971  PCLIP’CLIPP 

IF( PCLIP. LE.0. ) GO  TO  1000 
T0=TDL 
R=RJl 

CALL  TSVAUW.D.R.T0) 

1F( PCLIP. CE. PA( 1))  GO  TO  972 
TCC=-AA!  l)*AL0G!  (PCLIP- PC(  1)  )/PA(  D) 

PCL!P=PA( 1 ) *EXP! -TCC/AA! l) >+PC!  1)-PD( l)*SIN!P!*TCC/TTt 5) ) 

WRITE!  3. 38) PCLIP. TCC 

38  FORMAT!  1211  PCLIP! PSI) 3 ,FI3. 6. 10X.8R  TCLIP-.F13.6.4H  8EC) 

WRITE!  11,395  PC.L I P , TCC 

39  FORMAT! 28  //.  12H  PCLIP(PSl)3  .F13.6,  10X.3H  TCL1P" ,F13.B, 

'4H  SEC//) 

CO  TO  973 

972  TCC-O. 

PCLIP=0. 

973  TC-TO+TCC 
PASS* - 1 . 

WRITE!  11,42) 

42  FORMAT! 22H  REFL.  SHOT  PARAMETERS) 

IF1RD.LE.0.0)  GO  TO  63 
GO  TO  64 

1000  ETST=CEN*EMAX* 1 . E-8 
ENDFILE  It 
C 
C 

C OUTPUT  COMPUTATION  RESULTS 

DO  910  1= l.NFREQ 
XL= ! l-l) 

FREQ=  XL*DFREQ+FO 

ENRG=CKN*!RKFP!  I)*REFP!  D+XIMFP!  1)*XIMFP(  I> ) 

I F!  ENRC . LT.  ETST)  ENRG=  ETST 
XLEVL= 10. *AL0C10!  ENRG) 

910  WRITE!  1 1 .20)  FREti.REFP!  D.X1MFP!  1 ) . ENRC . XLEVL 
ENDFILE  It 
WRITE! 3, 21) 

3 FORMAT!  1 HI. 30H  FOURIER  SPECTRUM  OF  SHOT  DATA) 

6 FORMAT!/ 1 2H  YIELD! LBS) = ,F8.2/ l IK  DEPTH! FT) * . 

'F9.2UU  RANGE!  FT)  3 , F9 . 2) 

8 FORMAT! /23H  CONSTANTS  OF  TIME  SERIES/, 

'4511  IN  FORM  PRESSURE-TIME  CONSTANT  FOR  EACH  TERM/, 

'711  PERIOD. 8X.  2HPA, 8X, 9U! THETA  A)  8X.2HPB.6X,  911!  THETA  B)  . 

' 10X,  21U’C , 9X,  211PD) 

9 FORMAT!  IX, 2!  F8.5.3X)  , 111!  ,E9.4,  IH)  .OX. 

'F8.5.3X,  III!  , E9 . 4 . Ill)  .2F13.5) 

10  FORMAT!  911  INFINITE. Fl  1.3. 3X,  UK  ,E9. 4.  UD/llil) 

7 FORMAT!  1 9H  STARTING  FREQ!  IP/.)  » , F 10 . 2/ 

'1211  NO.  OF  FRF,U3  . 3X.  F 12 . 0/ 

'1911  FRRU  1 NTF.RV  AL!  H7.)  3 . F 1 0 . 2> 

12  FORMAT!  34H  FREQUENCY  SPECTRUM!  1»S  1 'HZ>  ENERGY  ENERGY  LEVEL/ 
'3X.4I1!  It2>  .6X,  1411REAL  1MAC.4X.23H! EHC/CM2/HZ)  < DB  RE  1.)) 

20  FORMAT!  F 50. 2.  IX.  2F.lt).  4 , K(2. 6 . F 14 . 2) 


173  21  FORMAT!  01 H DATA  HAS  BEEN  TRANSFERRED  TO  NT0!> 

176  40  FORMATt  lit  . 12HT»(SEC)=  .F13.8) 

177  GO  TO  15 

178  1010  STOP 

179  FND 

O ERRORS  COMPILATION  COMPLETE 


Table  A- 2 Con't 


t c 

2 COW  LEX  FUNCTION  U2(THlN,OM.El.T> 

3 C COMPUTES  COMPLEX  VALU  OF  2TUEXP( lOMELT)xT82+0M2) 

4 COMPLEX  K( . ARC 

5 

6 DEN'THIN*THlN+OM*OM 

7 U2»2.*TBlN*CEXP(El*OM*ELT)/DEN 

8 RETURN 

4 END 

• ERRORS  COMPILATION  COMPLETE 


1 C 

2 COMPLEX  FUNCTION  U3(0H.ELT,Dtm 

3 C COMPUTES  VALU  OF  Cl  EXPU  OM  ELT)  <£30*11  OM  DUT)*i  )ACI2-0K2) 

4 COMPLEX  El. ARC 

5 DATA  P l. El. 0/3. 14 1392684. <0. . t.> .3E-3/ 

6 ci*pi/uirr 

7 X*  C 1 -OM 

8 ADEL=>  ABS<  X) 

9 1 F ( ADEL . l.T.  O)  CO  TO  20 

10  ARC3  E l *QM 

U U3*C1*CEXP<  ARG*ELT1*(CEXP<  ARC* DUT>* Cl*OH>*X> 

12  CO  TO  30 

13  20  Xs  X*  OUT 

14  ARC-  ( Pf /<  Cl+OM) )*<  X*< 1 . -X»X/ 12. ' /2.  •*•£!»<  I.-X*X/6.)> 

13  U33CEXP<EI»0M*ELT>*ARG 

16  30  CONTINUE 

17  RETURN 

18  END 

0 ERRORS  COMPILATION  COMPLETE 
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Table  A- 2 Con't 


C 


C 


c 


c 


c 


< 


1 c 

2 SUBROUTINE  TSVALl W. D.R.TST) 

3 COMMON  T<3>.TT<4>.PA<3),PB(4>.PC<4>,PD<4>,AA15),AB(4) 

4 C PROGRAM  COMPUTES  VALUES  FOR  CONSTANTS  IN  TIME  SERIES 

3 PI«3. 141592*54 

* Z*D*33. 

7 VI*V**( I .<"3. ) 

8 XI=VI^Z»*<3./6.) 

9 X2*Vl/R 

!•  TT<  l)«4.34*XI 

11  TTt2)  *3.8**XI 

12  TTC  3) 32.48*X1 

13  TT(4)=2.3I*XI 

14  TU>»TST 

15  Tl 2) STT( I ) ♦Tt  I ) 

I*  T(3)*T1<2I*T(2) 

17  Tt  4>  *TT(  31  +T(  3) 

18  Tt  3)  =TTt  4)  ♦Tt  4> 

19  C TT  IS  THE  TIME  INTERVAL  BETWEEN  PULSE  PEAKS , T IS  THE  ELAPSED  TIME 

28  C TO  THE  PEAKS 

21  C THE  COMPONENT  PRESURE  AMPLITUDES  IN  EACH  TIKE  INTERVAL  ARE  PA.PB. 

22  C PC.ARL  FD 

23  POt I)=0.*(/**(2./3.i»»X2 

24  PCI  1 ) » . 483*PDt I ) 

25  PAl |)*2.88F.4*tX2**l. l3>-PCt 

26  P l - 3308 . «X2 

27  PBt l>aP!-PCl l> 

28  C THE  TINE  CONSTANT  FOR  EXPONENTIAL  TESItS  IS  AA 

29  AA(  l)®3.8E~3*W|*)C2**-.22 

3*  PITHl'PAt I ) *AA( I) 

31  PlMTtaPDt  I » *TTt  II 

32  Cla2.'P1-.483 

33  AA(2)MCU»PlKTl-PlTUI)/PBU> 

34  AB(D*AA<2> 

35  PA(2)'PBU> 

3*  PBt  2>  5 . 22#P  t-PCI  1> 

37  PC12>*PC<  l) 

311  C2=EXPl-.  I33*rri21/AA(2>  > 

39  PIH  2>  *2. 42*t PC« 2> ♦PAl 2) «C2* 

4#  C3S  2 . «PB<  2 1 ' P I -PCt  2> 

41  AAt3>«tC3«TTi2>'TAt2J»AAt2>  »-'P»t2> 

42  ABt  2* a AAt  3* 

43  PA13>*.22»1M 

44  PCt3»=60 

43  PBt3>*.IO*PI 

4*  PAI4MPIM3T 

47  PBt4M  ,03*Pl 

48  PAI35»S*IM4» 

49  P0t4>*0  8 

38  AA(  5' * \At  3> 

31  AAt  4 • * AAt  3 5 

32  ABt  3 > • AA<  4 > 

33  ABt  41 » AAt  3> 

34  Cl*  . U.*PI*IM«AA«3> 

S3  PDl3>»CI/'TTt3l 

3*  P&t  4'  * 8 . 3»C  I ■/'TTt  4' 

37  RETURN 

58  END 


( 
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Table  A- 3 OUTPUT  OF  SHOT  SPECTRUM 


YIELD!  LB8)«  1.69 

DEPTH!  FT)  s 69.99 
RANGE!  FT)  • 3M.M 


CONSTANTS  OF  TIME  SERIES 

IN  FORM  PRESSURK-TI HE  CONSTANT  FOR  EACH  TERM 


PERIOD 

PA 

t THETA  A) 

PB 

(THETA  B) 

PC 

. 12083 

40.4478* 

( .2370E-03) 

13. I 1230 

( .6933E-03) 

.2*834 

.083  14 

13. 1 1230 

» .6933E-03) 

2.67344 

( . I222K-02) 

.2*834 

.«*90& 

2.44374 

( . I222K-02) 

1 . 3UU08 

< . 1222E-02) 

0 . 00000 

. 0643 1 

1. 33808 

! . 1222E-02) 

.40143 

( . I222K-02) 

0.00000 

INFINITE 

.401 

( . I222E-02) 

PCLIPtPSI) 

= 20.097213 

TCLIP- 

.0001*132  SEC 

REFL.  SHOT  PARAMETERS 

CONSTANTS 

OF  TIME  SERIES 

IN  FORM  PRESSURE-TIME 

CONSTANT  FOR 

EACH  TERM 

PERIOD 

PA 

(THETA  A) 

PH 

* THETA  8) 

PC 

. 12083 

40.4478* 

( . 2370E - 03  * 

13. 1 1230 

( .6933K-03* 

.2*834 

. 083 14 

13. I 1230 

< . *9331.-03 1 

2.67344 

( . I222E-02) 

. 2*834 

.06405 

2.44374 

( . I222K-02V 

l . 33808 

« . I222E-02) 

0.00000 

.0*431 

1 . 33808 

( . 1222E-02) 

.40143 

( . 1222E-02) 

w • WV«  V 

INFINITE 

.401 

( . 1222K-02) 

1 


FREQUENCY 

SPECTRUM!  PS \s HZ) 

ENERGY  ENERGY 

I.EVF.L 

IKZ> 

REA1. 

1 MAO  ( 

ERC-'CM2-'ttZ>  t DU 

RE  1 . > 

0.00 

. I447E-02 

OO'  >F.  01 

. I29273K  00 

-888 

1.00 

. I439F.-02- 

76V3E-04 

. 12*54411  OO 

-8.98 

2.00 

. I392E-02- 

5 I64K-04 

. 1 14730E  00 

-9 . 22 

3.00 

. I246F.-02 

I325E-04 

. 103*8*1:  OO 

-9.84 

4.00 

. 1 034 K- 02 

6308E-04 

.6*238 IE- 01 

-11.79 

3.00 

. * I77E-93- 

38*41-04 

. 237628K-0 1 

- 16.24 

6 . OO 

. 802OK-O4 • 

470l.lF.-03 

. 1407432-01 

- ia.3t 

7.00 

- . 3* 14K-03- 

13 14K-02 

. 144 332 C 0« 

-8.23 

8.00 

. I343E-83- 

31 I7E-02 

.On  lo;l3K  0(» 

-2.21 

4.00 

.2064E-02- 

4083E-02 

. 1241 HOE  01 

5 . II 

10.00 

•4470K-U2 

3304V- 02 

. I4087OE  >1 

281 

II  . (Ml 

.3746K-02- 

48 IIIE- 03 

. 2*874  * } »| 

3 . 20 

12.00 

. 384* F- 02 

. 24*7E-t»2 

. 145*177  O* 

l . *3 

13.00 

- . 83  4 HE  03 

. 2847V- **2 

. 3434*3 V.  OO 

-2-63 

14.00 

-.27141-02- 

i*rm>«*2 

. O30I64K  OO 

-2.01 

13.00 

. 282*1.-03- 

. 4*mV-02 

. 1313801  Ol 

1 . 14 

16.00 

.271 TV -02- 

.33321-02 

. I2253*F  ot 

.88 

»7.o» 

. 27*41.-02- 

2784V  02 

,433 182 E OO 

.21 

: 0 . oo 

• 3744V. -02- 

. 29*37 V-  02 

. 1 39768V  Ol 

1 . *5 

14.00 

. 4447V-02- 

. 131 OF- 02 

. 1*47:* IK  Ol 

2.  17 

20 . t»£» 

.4203E-02 

. 1 307  V- 02 

. 1 I3277K  01 

. *2 

■»  1 OO 

. IU33K-02 

.2146V- 02 

. 442036 ¥ oo 

-3.08 

22  t»o 

-. 12331 -02 

. I243K-02 

. 14*R,  7 V.  oo 

-7 . 06 

23  OO 

- . 30*0!  -02- 

. 2300V-0? 

4 043*3 V oa 

- . 44 

24 . OO 

- 10*31.-02- 

. *.ncir-«2 

274453E  oj 

4.38 

23 . OO 

. 34ft3t>«2- 

.7*3::*.-o2 

. 437404E  Ol 

* . 60 

26.00 

7{»i«tt-02- 

.3878}:- 02 

. 4?*HHOK  o? 

6.78 

27 . oo 

.65H3E-02 

. ( 240V  - 07 

.27*4  lltt.  Ol 

4.44 

28 . ‘*0 

-2J48E-02 

. 2047F-O2 

. 556(1 2K.  oo 

-2.53 

24.  os 

.83441-03- 

.2494V-03 

.4*8*  I2F.-OI 

- 13.30 

30.00 

. 1 1321.-02- 

.3070K-04 

. 740943*. -O , 

- n . «2 

31.03 

- . I82HK-02- 

. 23437-03 

.204740K  OO 

-6.78 

:i2 . *h> 

- .32OIE-02- 

. S-«3*V-02 

. 233707 F.  oi 

8.  *9 

3.1 . oo 

. 144*)  1-02 

. 4200F-02 

. S36293H  0 1 

7.29 

34-00 

.70411-02- 

. 7383V- 02 

.6-*  2334V  Ol 

H.o« 

33 . OO 

. 4 1 4OK-02- 

. 2*06  V>  02 

. 3* VP24E  Ol 

7.31 

36  <U1 

-7443V- 02 

23* IF  02 

. 0:7202 V.  PI 

3 38 

PD 

.6*58* 
. *4939 
. I IB99 
.0*387 


PD 

.**38* 
. *4939 
. 1 1849 
.06387 
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A.  3 The  Pilter  Time  Series  and  SHOTTIME 

As  described  in  Section  3 of  this  report,  the  time  series 
must  be  anti-alias  filtered  to  be  processed  using  digital  tech- 
niques. An  analytic  filter  can  be  applied  to  an  analytic  time 
function  through  a multiplication  of  the  complex  spectrum  of  the 
time  series  by  the  complex  filter  transfer  function,  and  trans- 
forming back  to  the  time  domain.  The  expressions  used  for  this  are; 


SQ(w)  = H(w)  Si(w) 


CO 

pc(t)  = '-j--  J H (co)  (oj)  e"’iaJtdu) 


(A. 2) 


(A. 3) 


A four-pole  Butterworth  filter  was  chosen  as  the  anti-alias 
filter  becav.se; 

a)  It  has  a flat  frequency  response  below  0,5  x f (cut- 
off  frequency) , thus  tr.inizing  the  need  for  frequency 
response  correction. 

b)  It  is  physically  realizable,  which  eliminates  the 
need  for  arbitrary  decisions  concerning  non-causal 
response. 

c)  Its  properties  are  analytically  well  known  and 
reasonably  tractable. 

d)  Its  physical  realization  is  readily  available 
commercially  and  is  frequently  used  in  this 
application. 
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The  squared  magnitude  of  the  frequency  response  is; 


2 8 
H(co)  | = 1/(1  + (^)  ) 


The  complex  frequency  response  is; 

4 i (w/m  ) 3 (w/w  )2 

(J±_\  4.  C C 

mc  sin  (n/8) 


H(ui)  = 1/ 


i (w/w  ) 


2 sin  (n/8) 


+ 1 


sin  (tt/8) 

This  function  has  four  poles  at  the  frequencies  given  by; 

•/.  . -e1,,/8,  -o-3"/a'  -o’5*/*.  -6-7’/8 


The  five  terms  of  H (co ) times  the  twelve  non-zero  terms 
of  S . ( m ) , when  integrated  using  standard  residue  techniques, 
result  in  a sixty  terra  expression  for  pQ(t)  which  may  be  inter- 
preted as  "steady-state"  response  to  the  applied  signal  plus 
transient  solutions  (i.e.  "ringing"  at  the  filter  poles)  assoc- 
iated with  the  starting  or  stopping  of  each  signal  component. 

As  a result,  the  significant  terms  in  each  epoch  may  be  easily 
selected  on  causal  or  significance  qrounds. 

The  interactive  Fortran  program  which  was  used  to  evaluate 
the  numerical  magnitude  of  each  term  is  called  SHOTTIME.  A 
listing  of  the  source  program  is  given  in  Table  A-4,  while 
a sample  output  listing  is  given  in  Table  A-5. 


Table  A- 4 LISTING  OF  SHOT  TIME 


i 

•> 

a 

4 

5 

6 

7 

8 
9 

10 
1 1 
12 
ia 

14 

15 

16 

17 

18 

19 

20 
21 
22 
25 

24 

25 

26 

27 

28 

29 

30 
2! 
52 
88 

84 

85 

86 
87 

au 

39 

40 

41 

42 

48 

44 

45 

46 

47 

48 
44 
50 

31 

52 

53 
34 
53 

56 

57 

58 


C NAME  SHOTTIME  3 ********  RC= 32000  ************ 

C COMPUTES  VALUES  OF  ANALYTIC  8H0T  TIME  SERIES  AT  UNIFORM  TIME 
C INCREMENTS  DDT  STARTING  FROM  TO.  THE  ARRIVAL  TIME  OF  THE  SHOCK 

C WAVE  IS  TIS.  MODEL  IS  MODIFIED  WESTON  (HALF  SINE  PLUS  CONST. 

C REPLACES  NEC AT I VE  CONSTANT  OF  WESTON  MODEL.) 

C INCLUDES  DELAYED  AND  ATTENUATED  SURFACE  REFLECTION  WHEN  NR  « 1, 

C AND  TRUNCATES  REFLECTION  AT  -CL1PC  WHEN  NR  > i.  DELAY  AND  ATTEN 

C ARE  INPUT  VALUES, OUTPUT  ARE  FILES  OF  500  VALUES  ON  11. FOR  TABLES t 

C FILES  OF  500  VALUES  + LAST  VALUE  REPEATED  FOR  PLOTS  ON  12 
C OUTPUT  ON  12  IN  BINARY  FORMAT 

C WHERE  26  IS  USED  AS  DUMMY  VARIABLE  NUMBER 

COMMON  T< 5) ,TT( 4) , PA<  5) ,PB(4) ,PC(4>  ,PD(4) ,AAi5) ,AB<4) 

COMMON  P 1 ( 500) , T9< 1) 

DIMENSION  N( 5) , TL( 20) 

P I a 3 . 141592654 
C INPUT  SHOT  DATA 
100  WRITE  (3, 1) 

wc-o.o 

READ(  2, 2) W, D. R1 . NR 

1 FORMAT! 37H  GIVE  YIELD(LBS) , DEPTH!  FT) .RANGE!  FT) , 

*19H  REFLECTION  NUMBER  /2l  H SEPARATED  BY  COMMAS:) 

2 FORMAT! 3F 13. 6. 12) 

I RCNT=  0 

IF!  NR.  EQ..  0)  GO  TO  300 
WRITE(3,3) 

READ!  2.2)  TDF.L , ATTEN 

3 FORMAT!  33H  GIVE  DELAY!  MSEC) , ATTEN  RE  UNITY:) 

T2S-TDEL/1000. 

IFINR.LT.2)  GO  TO  300 
WRITE!  3, 4) 

READ! 2. 2) CL  I PC 

4 FORMAT!  4411  GIVE  CLIP  CONSTANT!  PS  I ) FOR  NEGATIVE  PHASE:) 

300  WRITE!  3, 5) 

READ(2,2)T0,DT.XN 

3 FORMAT! 56 H GIVE  START  TIME! MSEC) .SAMPLE  TIME  INTERVAL! MSEC) .AND  NO 
*22H  OF  TIME  SERIES  VALUES/2 lH  SEPARATED  BY  COMMAS:) 

WRITE!  3.36) 

80  FORMAT!  44H  MOUNT  DATA  TAPE  ON  ITT© I. CURVE  TAPE  ON  HT02.) 

NRV.\L= NR 
Tt-VO/  1000. 

DDT=DT  1000. 

N70T=iXN*.5> 

T1S=0. 0 
UMAX” HTOT/50© 

I F(  ulMAX*500- NTOT) . NE . 0) JMAX= JMAX+ l 
DO  310  1= 1 . JMAX 

310  TL(  I ) - FLOAT! ! I- 1 > *300) *DDT+Tt 
PAUSE  l 1 

DO  1000  JT-  1 ,.1MAX 
AMP- I . 

NR3 NRVAL 
Tl-TL(JT) 

DO  330  l3 1,300 

ph  n = ». 

350  T91  D-O. 

NTOT 1=500 

l F!  .IT.  EQ . J MAX)  NTOT  1 = NTOT- ! J MAX-  l ) *500 
CALL  TSVAL(W.O.RI.TIS) 
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see 


660 


roe 


CALL  PTS! Tl , DDT, NT0Y1 , AMP. WC> 

OUTPUT  HEAD l NG 
WR1TF.(3,25)NT0T1 
IF! JT.GT.  l)GO  IX)  760 
WRITE!  11,6) 

WRITE!  II. 7) V.D. HI 
FG*WC'(2.*Pl> 

WRITE!  lt,8)r0,DT,XN,FC 
WRITE!  1 t ,2D 
CONTINUE 
WRITE! 1 1.22) 

DO  666  1-1,4 
L=  1+1 

WRITE!  1 1,0)  (Tl  1)  ,PA!  I)  . AA<  1)  ,PB(  l)  ,T(L)  ,AR(  I)  .PC!  I)  ,PD(I)  ,TT!  1)  ) 
CONTINUE 

WRITE! 11, 10)T(S) ,PA!5) ,AA!5) 

IF! NR.GT.0)  GO  TO  HOO 
IF!  JT.CT.  I MM)  IX)  060 
WRITE!  11,11) 

ENDFILE  11 

6 FORMAT!  1111,2511  TIME  SERIES  OF  SHOT  DATA) 

7 F0RHATt//l2H  YIELD! LbS) « , F8. 2/1 1H  DEPTH!  FT) * , F9 . 2' 1 IHRANGE! FT) • , 
*F9.2) 

« FORMAT!  / /U7H  START  TIME!  RE  SHOCK  DIRECT  ARRIVAL)3  , F0.5.6H!  ttSEC)  / 

• 17H  SAMPLE  INTERVAL®  , F9 . 5 ,611!  ffcSEC)  '/  1511  NO  OF  SAMPLES® .F9.0, 

*10X.  1311  FILTER  CUTOFF* , F 13. 4, 4H  RZ.  ) 

21  FORMAT! //51H  CONSTANTS  OF  TIME  SERIES.  TIME  SERIES  GENERAL  FORK/ 
*6211  P(Tl)*PA*EXP!  I T-TI ) MU)  + PH* EXP!  1 TI-T1)  /AB+PC-PD*S IN!  Tl-T)/TD  ) 

22  FORMAT!  ^/'  SX.ellTI  SEC)  , 5X,  7HPA1  PS  l ) . 5X.  7HAA!  SEC)  ,5X,7HPB!PS1)  ,5X. 
*711T1 1 SEC) , 5X, 7HAIH SEC)  . 5X.7HPC!  PSl)  , 5X.7HPD!  PSD  ,5X,7HT1'!8EC)) 

4 FORMAT!  F12.8.F12.6.F12. 10, F12.6 , P 12. 8, F 12. 10.2F12.6,F12.8) 

10  FORMAT!  F 12 . 8 , F 12  > 6 , F 12 . l0.4X,3H0.0,6X.8HlMFlNlT]i:  .3X.3H0.0. 13X, 

*3110 . 0 . 9X.  3H6 . 0 . 7X.  3110 . 0 - V ) 

11  FORMAT!  1HI.//50X,  U1HT1ME  SERIES  VALUES// 

*511  NO.  , 3 1 26 HT 1 ME!  MSEC ) PRES!  PSD  ,3X)  //) 

900  CONTINUE 

OUTPUT  TIME  SERIES  VALUES  IN  PACES  OF  160  LINES  WITH 
5 SETS  OF  VALUES  PER  LINE,  FOR  HARD  OOPY 

FOR  PLOT,  300  VAI.UES  FOR  1ST  PACE.  301  OR  LESS  FOR  CONTINUING 
PACES  WITH  LAST  VALUE  REPEATED  AT  TOP  OF  NEW  PACE,  IN  BINARY 
OUTPUT  DATA  FILES 
Nt>  ! JT-  1 ) *306 
DO  020  KM.NT0T1.5 
Nt  ' K+NG 
N2  * K+4 
IK)  9 tO  DK.N2 
T9l  I ) * 1000. *T0!  D 

WRITE!  I 1.23)  Nt,!T9<  l)  .Pit  l)  , DK,N2> 

CONTINUE 

FORMAT!  (5 . 5t  F 10. 4, FlO.6 , 3X» ) 

CO  IX)  1606 


910 


83!) 


BRANCH  W HERE  FOR  NR.CT.0.  COMPUTE  AND  ADD  IN  EFFECT  OF 
SURFACE  REFLECTION  IF  NR® 1 . TUUNCVTKS  VALUES  AT  -OH PC  WHEN  NR. 
COMPUTE  NEW  TSVAl.  AND  PRESSURE  TIME  SERIES.  IS  VALUES  ARE 
SUBTRACTED  Fll!)M  1*9. 
tUt-  Rl  - ATTKN 

CALL  ISVALi  W,  D.RR.T2S) 


> l 
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IFICLIPG.GT.PA!  1)  )NR*  1 
1 F!  NR. LT. 2) GO  TO  801 
AMP=-CL1PG 
GO  TO  802 

801  AMP*-2.0*PA(  i) 

802  CALL  PTS! Tl , DDT, NT0T1 ( AMP, WC) 

WRITE!3,735> 

7S5  FORMAT!  11 H OK  AT  * 5) 

1F( JT.CT. 1)G0  TO  900 
I Ft  NR. LT. 2) GO  TO  840 
WRITE! I 1 ,311  AMP 

81  FORMAT!  34H  THE  NEGATIVE  PHASE  OF  'IUE  DELAYED  PLUS  REFLECT  WAVE  , 

* 18UHAS  BEEN  TRUNCATED/ 13H  TO  THE  VALUE, F 1 1 . 6 ,3H( PS1 ) ///) 

840  CONTINUE 

1F< JT.GT. 1 )G0  TO  900 
ATDB320. *ALOG10!  ATTEN) 

SATDB- 1 , 13*ATDB 

WRITE! 11, 32) RR.TDEL, ATTEN ,ATDB, SATDB 
32  FORMAT! 36B  THE  SURFACE  REFLECTED  WAVE  AMPLITUDES  WERE  COMPUTED  FRO 
*24BM  A NEW  SET  OF  CONSTANTS/ 

*57B  THE  CONSTANTS  ARE  BASED  ON  A RANGE  FOR  THE  REFLECTION  OF, FI 0.2 
*,4H!FT)/15H  AND  A DELAY  OF,F!0.6,6H!MSEC)// 

#38H  THIS  CORRESPONDS  TO  AN  ATTENUATION  OF, F9.7, 1H! .F6.2.4H  DB,/ 
#51H  EXCEPT  FOR  THE  SHOCK  WAVE  WHICH  18  A'lTENUATED  BY.F6.2.4B  DB) 
#//58H  THE  NEW  TIME  SERIES  CONSTANTS  FOR  THE  REFLECTED  WAVE  ARE://) 
NR=0 

GO  TO  300 

C RETURNED  TO  PRINT  TIME  CONSTANTS  AND  TIME  SERIES  VALUES 

1000  CONTINUE 
ENDFILE  1 1 
ENDFILE  26 
WRITE!  3, 100!) 

1001  FORMAT!  2 1U  NEW  RUN? 1 , a YES 1 0 , * NO) 

READ!  2.23) l 
IF!  1.  Eft.  I)  CO  TO  100 
STOP 
END 

ERRORS  COMPILATION  COMPLETE 


f 


« 
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2 

3 C 

4 

5 
4 
7 
0 
4 

10 

11 

12 

13 

14 

13 

14 
17 

10  C 
14  C 
20  C 
2t  C 
22 

23 

24 

23 
26 

27  G 
20 

24 

30 

31 

32 

33 

34 

33 

36 

37 

30 

34 

40 

41 

42 

43 

44 

43 

46 

47 
40 

44 
34 

31 


Table  A- 4 Con*t. 


SUBROUTINE  TSVALl V.D.R.TST) 

COMMON  T(B) ,TT<4> ,PA<5> ,PB(4) ,PC(4) ,PD<4) ,AA(5) ,AB<4) 

PROGRAM  COMPUTES  VALUES  FOR  CONSTANTS  IN  TIME  SERIES 
Pl*3, 141342634 
Z»  D+33 . 

Wl*  W**< 1./3.) 

Xl»Wl/7**<5./6.) 

X2-V1/R 

TT(  1 ) =4. 34#Xl 

TT<2)a3.06*XI 

TT<  3) =2. 48*Xl 

TT(  4)  =2.31#X1 

Tl  1 ) *TST 

r(2i*TT(  1)+Tl  1) 

T(  3)  = TT(  2)  +T<  2) 

T(  4>  =TT<  3)  ♦Tl  3> 

TlB)=TT(4)+T(4) 

TT  IS  THE  TIME  INTERVAL  BETWEEN  PULSE  PEAKS, T IS  THE  ELAPSED  TIME 
TO  THE  PEAKS 

THE  COMPONENT  FRESURE  AMPLITUDES  IN  EACH  TIME  INTERVAL  ARE  PA,PB, 
PC. AND  PB 

PD< 1>«B.*<Z**<2./3.)>*X2 
PC( l>  = 403*PD( 1) 

PA( t)«2.0BK4*<X2**l. I3)-PC(  t) 

P 1=3300. *X2 
PB< D-P1-PC1 1) 

THE  TIKE  CONSTANT  FOR  EXPONENTIAL  TERMS  IS  AA 
AA1  0=3. 8E-5* V 1 *X2**- . 22 
PlTHfFAl  l)»AA(  1) 

PlKri  = PD<  1)»TT<  1) 

Ct=2. /PI- . 403 

AA(2)  = (Ct*Plffri-PlTHl)/PB(  1) 

AB(  1 ) = AA(  2) 

PA(2>=PBt l) 

PB1 2) « .22*P l-PCl t) 

PCI  2) = PCI  1) 

C2SKXP< - . 135*TTv  2)  ^ AAl  2) ) 

PD( 2) = 2 . 42*1  PCI  2) ♦PAl 2) *C2> 

C3*2.*PD(2)/Pl-PC(2) 

AAl  3) « 1 C3*rr<  2)  -PA(  2)  *AA(  2)  > /PB1  2> 

ABl  2)  = AAl  3) 

PA13)= .22*P1 
PCI  3) =0.0 
PB<  3)  = . 10»Pl 
PA1 4) * PBt  3) 

PB( 4) * . 03*P 1 
PAl 3) *PUi  4) 

PC(4)°O.0 
AA( 3) 0 AA( 3) 

AA<  4) * AAl  3) 

ABl  3' * AAt  4) 


32  ABl  4) » AAi 5> 

33  Cl" . 16*Pl*Pt*AA13) 

34  PD(  3)  "Cl  ''TIT 3) 

83  PDl4)=0.3»Cl/TT(4) 

36  RETURN 

37  END 

0 ERRORS  COMPILATION  COMPLETE 
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Table  A- 4 Con't. 


1 

2 

3 

4 

5 
4 

7 

8 
9 

10 
1 1 
12 
r.i 
14 

is 

16 

17 

18 

19 

20 

21 

<)<i 

2~\ 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 
4« 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 
3 I 

52 

33 

34 

53 

56 

57 
5U 


SUBROlTf  1 WE  PTS!  T0, DDT, NTOT, AMP , VC) 

COMMON  T(S) , TTC4) , PA(5) ,PB<4> ,PC(4) ,PD!4> ,AA(5) , AB(4) ,Pl<5©0) 
COMMON  T9( 1) 

C COMPUTES  PRESSURE  TIME  SERIES  VALUES  FROM  T0  III  STEPS  DDT  FOR  NTOT 
C STEPS,  TIME  AND  PRESSURE  VALUES  RETURNED  IN  T9  AND  PI 

C VERSION  2-  SIGNAL  LOW  PASS  FILTERED  THROUGH  A 4P0LF  BUTTERLORTH 

C (48DB/0CTAVE  SLOPE)  AND  -3DB  POINT  AT  WC  FILTER  IS  IMPLEMENTED 

C ANALYTICALLY  ON  THE  ORIGINAL  ANALYTIC  TirtE  SERIES  EXPRESSION 

COMPLEX  0M(2) ,HR(3,2) , El , A1 , A2,HRN< 2) .ARG.CP.Bl ,B2,HI,H2,H3 
COMPLEX  EX.H,PIH,WR(2) ,CPP 
DIMENSION  R( 4) , A( 4) 

DATA  El.  Al,  A2A0. .!.),(  .9238795, . 38268343) ,(  .382683431.  .9238795)/ 
DATA  B1.B2/I-. 191341 .- .46 19397) . ( l . 1 15221 . .46 19398)/, PI/3. 141593/ 
OM.0E-8 

IF( WC.NE.0. ) GO  TO  8 
WR!TE(3,2) 

2 FORMAT! 3511  GIVE  FILTER  CUTOFF  FREO. ( HZ, F12.4) ) 

READ! 2.3) WC 

3 FORMAT!  F 12 . 4) 

WC-2. *P 1*WC 

C FIND  CL I PlINC  PARAMETERS 

8  IF!  AMP.GT.0.)CO  TO  9 
PCLIP=-AMP 

IF!  PCHP.  GT.  PA!  l))GO  TO  6 
TCC=-AA!  l)*ALOG( ! PCLIP-PC!  1))/PA!  D) 

PCL1P=  PA!  1 ) *EXP!  -TCC/AA!  1))+PC! 1)-PD! 1) *SIN! PI»TCC/TT! 1)) 
TC=TCC+Tt  1) 

6 AMP1- 1 

9  OM!  l)=WC*-t.*Al 
0M!2)=WC*-1.*A2 
HRN! I ) =WC*BI 
HRN!2)*WC*B2 
DO  10  1=1.4 
J«l4| 

R!  D=PI/TT!  I) 

A!  I ) = l . / AA!  J ) 

10  CONTINUE 
DO  20  N= 1 .2 

HR!  1 ,N)=HRN!N)*(PA!  l)/!OM!N)+El/AA!  I)) 

♦♦PC!  1 ) /ON!  N)-EI*PINt  PD(  1)  ,OM!N)  ,R!  1))) 

HR!2.N)=HRNCN»*Et*-l  .«!PIN!PD!  I)  .ORtN)  . R<  l>) 

♦♦PIN!  PR! 2)  , OM!  N)  , R!  2)  >+2.*PA!2)*At  l >/!  OH!  N>  *OM!  N) +A!  l)*A!  1))> 

HR!  3 . N > = HRN ! N ) *!  P I M!  PD!  2)  ,OM(N>  ,IU2> ) 

♦+P 1 Ml  PD!  3)  , ON<  N)  , Ri  3) ) ♦2.*PB<  2)  *A(  2»/l  ON!  N)  *OM!  N)  + A1 2)  *A!  2)  ) 
♦♦El*PC12>/!<m<  N ) + E I / AA! 3 > )>*KI*-1. 

UR!4,N»=HUMN)*EI*-l  .*(IMN!  PDi3)  ,OM(N)  ,R!3)  > 

*+P 1 M!  PI)!  4)  . DM.  N>  . R(  4)  ) +2. *PAi  4)  #A<  3)  /!  OM!  N)  »OM<  N)  ♦A!  3)  *A! 3) ) ) 

HR!  3.N)  * HRN!  N) »!  P IN!  PDl  4)  .OMlN)  , RC  4 ' ' 

♦♦2. *PAC 3 t*A! 4) /) OM! N) *OM!  N) ♦ At  4)* V 4) ) > *El»- 1 . 

29  CONTINUE 

DO  30  1= 1 .3 
DO  30  N*  I .2 
HR!  I . N > * 2 . * t UR! 1 ,N) ) 

30  CONTINUE 

OM!  I > = E I *0M<  I ) 

ON! 2t  ElvOMtZ) 

4 FORMAT!  4E 12. 4) 

C TEST  LOOP  PARAMETER 


i 

i 
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Table  A- 4 Con't. 


04 

IF!  DDT.  GT.  9.)  CO  TO  09 

6# 

VR1TE13, l) 

61 

1 

FORMAT! 0 1 H FATAL  INCREMENT  ERROR  lit  8UB  PTB, ABORT  RUM  AMD  FIX) 

62 

701 

FORMAT!  1 1 H OK  AT  * 1) 

63 

702 

FORMAT  ( 1 IB  OK  AT  e 2) 

64 

753 

FORMAT! 1 IB  OK  AT  • 3) 

60 

754 

FORMAT! 1 IB  OK  AT  * 4) 

66 

GO  TO  1999 

67 

C 

INITIALIZE  VALDES 

«B 

09 

J*1 

69 

1«1 

79 

TI»T0 

71 

IFtTI.GE.TlO)) J*8 

72 

TMI*T(  1) 

73 

THAXH* FLOAT!  NTOT) «DDT+T9 

74 

99 

L»J+1 

75 

1F!L.E0.6)G0  TO  499 

76 

TMAX«T!L) 

77 

IF!  TMAX.  CT.  TMAXM)  TMAX*  TMAXH 

78 

1 F( T I . LE . TKAX) GO  TO  199 

74 

J»JM 

86 

CO  TO  99 

at 

199 

COMTINUE 

82 

C 

CHECK  FUNCTION  IN  RANCE  OR  COMPUTE  SPECIAL  VALUES 

83 

IFlTt.CE.TMDGO  TO  299 

84 

Pl(  D-9.+P11  I) 

85 

T9!  1)»TI 

86 

TI»DDT*FLOAT!  D+TO 

87 

I*  !♦! 

88 

CO  TO  199 

89 

299 

IK(Tl .GE.Tt 8) )GO  TO  499 

99 

C 

COMPOTE  FUNCTION  VALUES 

91 

209 

ARC*E1*-1./AA!J) 

92 

Hl-H( ARG,VC) 

93 

ARC»E1/AA!L) 

94 

H2-H« ARC, VC) 

90 

ARG"Pl/'TT!  J) 

96 

H3*H( ARC, VC) 

97 

ARG-1PI/1T!  J))*El 

98 

399 

V»(T!  JJ-TD/AA!  J) 

99 

V«EXP!V)»PA!  J) 

199 

B»!T1-T!  L))/AA!U 

191 

B*EXP!  B) *PB<  J) 

192 

391 

l?(Tl.GT.T!2))GO  TO  329 

193 

If!  AMP.CT.9. )GO  TO  329 

194 

IF! PCL1P. GT.PAt 1)>G0  TO  329 

190 

VR( l)«-Al»VC 

196 

VR(2)»-A2*VC 

1*7 

IF1TI.CT.TC)  GO  TO  319 

198 

CP»HRNt  t)*CEXPtEI«(Tt  l)-TI)«VRt  |))/VR(  p 

199 

CF»BIW!2l*CEXPtEl*tTt  P-TI>*VR!2>  >/VR(2)4CP 

119 

PI!  1 ) *PCLIP*<  l.42.9*REAL!CP))*AMP*PH  l) 

111 

CO  TO  349 

112 

319 

AARG*  P 1 •TCC/TTl  1) 

113 

CP*V*H14B»R2»PD!  J)*A1MAC!  H3«CEXPt  !T(  J)-TI)*ARG)  >*PCt  1) 

114 

DO  311  1 J-  1 .3 

110 

CPP»CEXP<  El*Vft(  U)*!T!  1)-T1)>*HRN!  IJ>*2.9 

114 

CP«CP4CPF»! PA!  li/tWlt  1J)*EI/AA!  imPC!  P«CEXPt £t»TCC«VRt  1J)>/ 
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Table  A- 4 Con't. 


117 

#WR<  1J)> 

na 

CPP=CPP*CEXF<E1*VR< lJ)*TGC>/iWR(  tJ)**2-(Pt/TT<  1))*42)«PD( 

1 1'» 

CP=CP-CPP*<E1*PI*C0S<AARC)/TT<  1>4WII(  IJ)*SIN< AARC) > 

120 

3 : i 

CONTINUE 

121 

312 

IM( |)-P1( U +AHP*REAL< CP) 

122 

CO  TO  340 

122 

320 

CP*  V*Hl+B*H2+PD<  J) *AIMAG(  H3*CEXP(  ( T(  J) -TI) *ARG) ) 

124 

*+RE4L<EX<Tl,0N(  11  ,T(  J)  >*HR<  J.  1>+EX(T1 ,0H(2)  ,T<  J)  )*HR(  J.2> 

125 

Pl(  1>  = (REAUCP>+PC<  J))*AMP4P1(  1) 

126 

340 

VH  1 ) *TI 

127 

T1=DDT*FL0AT<  D+TO 

120 

1=1+1 

12‘) 

I F(  T 1 . LT.THAX1GO  TO  300 

toe 

1F(  l.GT.NTOT)GO  TO  1000 

131 

J=J+1 

132 

L=J+I 

133 

TMAX=T(  L> 

134 

1 F(  TMAX.  CT . TWAXM)  TTtAX*  TNAXH 

133 

IF( J.l.T.S'GO  TO  280 

136 

400 

1F( 1 . CT. n IOT)  O'  TO  1000 

137 

ARG=El*-l.'-riAC&) 

13U 

HI*PA(5)  *m  ARC.  VC) 

130 

430 

DtF=T(5)-Tl 

140 

CP=  II 1 *£XP(  D I F*  A(  4>)  +RF.AL<  CEXP(  DIF*OM(  1 » ) *HR(  8.1) 

141 

3+CEXP ( 0 1 F*OM<  2>  1 «HR(  8,2)1 

142 

Pl<  ! > = REA1.<CP>*ANP+P1(  1) 

143 

T4(  I ) = T I 

144 

Tl  = DDT»FLOAT(  M+TO 

14S 

1*1+1 

146 

500 

IF< 1.CT.NT0T)  CO  TO  1000 

147 

IF<  ABS(  REAl.(CP)  ) .GT.ft)  GO  TO  430 

140 

N 1*  1 

140 

DO  600  l*Nl.NTOT 

130 

PH  l>*0.+Pl(  l* 

131 

T4(  l)=DDT* FLOAT!  D+TO 

132 

600 

CONTINUE 

133 

1O00 

1 F(  ANP . L T. 0 . ) ANP* -PCLIP 

104 

RETURN 

135 

END 

0 

ERRORS 

COMPILATION  COMPLETE 

r 


i 


[ 


! 

\ < 


i 
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Table  A- 4 Con*t. 


t COMPLEX  FUNCTION  H(A,W) 

2 COMPLEX  EI.A.C 

0 E!=<«.,1.) 

4 C=W 

3 D*SIN<3. 141993/8.) 

6 Hs  1 . ✓(  C*C*C*C-C*C/'<  2.  *D*D)  * 1 . +EI*( C*< C*C- 1 . ) ) /D) 

7 RETURN 

0 END 

» ERRORS  COMPILATION  COMPLETE 


1 COMPLEX  FUNCTION  PIMCP.O.T) 

2 COMPLEX  0 

3 PlM=P*T/(0*0-T#T) 

4 RETURN 

5 END 

0 ERRORS  COMPILATION  COMPLETE 


1 COMPLEX  FUNCTION  EX(T.O.A) 

2 COMPLEX  O 

3 EX=CEXPUA-T)*0> 

4 RETURN 

5 END 

0 ERRORS  COMPILATION  COMPLETE 
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Table  A-5  OUTPUT  OF  SHOT  TIME 
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APPENDIX  B 

ERRORS  ASSOCIATED  WITH 
REALIZABLE  ANTI-ALIASING  FILTERS 
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APPENDIX  B 


Errors  Associated  with  Realizable  Anti-aliasing  Filters 

Equations  3.1  and  3.3  provide  a basis  for  estimating  the  error 
due  to  small  amounts  of  aliasing  for  several  important  cases 
that  will  be  discussed  here.  will  consider  two  cases,  that 
where  the  energy  density  spectrum  is  constant  vs  frequency 
(corresponding  to  that  of  an  impulse  of  infinitesimal  duration) , 
and  that  where  the  spectrum  rolls  off  smoothly  at  6 dB/octave 
above  some  frequency  (corresponding  to  an  impulse  with  an  expo- 
nential decay).  For  both  cases  we  will  employ  an  ideal  n-pole 
Butterworth  anti-aliasing  filter  as  representative  of  a real- 
izable filter.  The  first  case  represents  a conservative  approach, 
which  is  applicable  for  an  arbitrary  shot  spectrum,  while  the 
second  case  may  bo  used  when  more  information  is  available  about 
the  spectrum. 

For  the  case  of  weak  aliasing,  we  can  treat  the  effect 
of  aliasing  as  noise  that  contributes  to  the  computed  spectrum, 
and  compute  an  equivalent  signal  to  noise  ratio. 

Prom  Equation  3.2  we  have  in  terns  of  energy, 

EMf)  - V 2P*  ( f - 2qf  ) 

“ * m 

qa-* 

where  f = 1/At 

tQ 

For  an  n-pole  Butterworth  filter  with  cutoff  frequency  f c , 


we  have 


H2(f) 


1 ♦ (f/fc) 


Tn 


Thus#  for  Case  1 (noting  that  the  energy  is  the  sum  of  equal 

2 2 

contributions  from  P (f)  and  P (-f)). 


Ex(f)  = 2PJ  (f ) = 


2C 


1 ♦ (f/fc) 


nr 


, f > 0 


(B.l) 


and  for  Case  2 (with  roll  off  at  f^) 


E2(f)  » 2P2(f)  = 


(B.  2) 


2C  

(i  + u/f  )2)  u + tf/f  )2n) 

r c 


,f  > 0 


(B,  2) 


CASE  1 - FLAT  SPECTRUM 

Substituting  B.l  in  3.4  anu  simplifying  yields: 

1 ♦ 


E,  (f)  = E.  (f)  + 2C  l 
1 q=l 


2n 


1 


n 


We  note  that  f _<  f and  make  the  reasonable  assumption  that 
f ^ f^.  Iqnoring  1 in  the  denominator  of  each  series 
(since  the  terms  in  parenthesis  are  greater  than  1),  we  have: 


a- 2 


E1  (f)  = E1(f)  + 


n > E1 


q=l  (a+q) 


(a-q) 


where  a = 


By  our  assumptions,  0 < a £ -j—,  so  we  can  replace  the  series 
by  the  following  upper  bounds: 

l — - — 2^  i l -?r=  5( 2n) 

q=l  (a+q) 2n  q=l  q2n 

00  00 

l ” — Yn  - £ 2n 

q=l  (a-q)2  q=l  (q-^)2 


* + 


+ l = 22n  + Stfn) 


q=l  q 


where  £ (n)  is  the  Rieman  Zera  Function,  for  integer  r..  Some 
values  of  this  function  are  given  in  Table  B-l. 

The  upper  bound  estimate  for  the  energy  spectrum  with 
aliasing  is, 

e!  (f)  = £.  (f)  + r--,  • (22n  + (2n) ) > E*  ' > e! 

1 1 / 2fn  V 1 1 


The  signal  to  alias-noise  ratio  is. 


E (f)  /2fN\ 

SN,  (n, f ) < — rr? = r3  /( 

A ~ E:  (f\-E}(f)  \rc  / 


(2*“+2C{2n)) (l+(f/f  ) ") 


M.  Abramowitz  and  I.  A.  Stegun,  Handbook  of  Mathematical 
Functions,  National  Bureau  of  Standards,  Applied  Mathematical 
Series,  Chapter  23.2  and  Table  23.3  (June  1964) 


Table  B-l 


Values  of  Riemann  Zeta  Function 


n 

Un) 

2 

1.64493 

4 

1.08232 

6 

1.01734 

8 

1.00407 

10 

1.00099 

12 

1.U0024 

(B.  3) 


SN1(nf)<(fN/fc) 


for  f,  f < f 
c — m 


2n/2 


Equation  B.3  may  be  converted  to  decibels  by  taking  10  log^ 
of  the  right  hand, side.  Table  B-2  provides  an  evaluation 
of  Equation  B. 3 for  various  values  of  f /f^  at  the  frequency 

f - V 


CASE  2 SPECTRUM  ROLLOFF  AT  f . For  this  case: 

E’(f)  = E (f)+2C  l 5— 

q-1  (l+((f+2qfN)/fr)^) (l+((f+2qfN)/fc)Jn) 


(l+((f+2qfN)/fr)2)  (l+{(f-2qfN)/fc)in) 

with  the  assumption  that  fr  2f^,  in  addition  to  the  assump- 
tions of  Case  1.  We  can  ultimately  reach  the  simplified 
expression  for  an,  upper  bound  to  the  aliased  energy  spectrum 


Ej  ({)  - Ej(f>  + -n— JTTT 

(-jS)  <-A 

r c 


with  signal  to  noise  ratio , 


(22n*2+2U2n+2)) 


2f  2 2f  2 
SU.4  a (~)  / (2 


2n+2 


+ 2£(2n+2) (1+f/f  )2) (1+f/f  )2n) 

X v 


f , 2 £ 2n 

^ ijr'-)  (»-“)  /4 

r c 

Study  of  this  expression  indicates  that  the  presence  of 
roll  off  of  6 dB/octave  commencing  at  a frequency  £ is  roughly 


Table  B-2 


SIGNAL  TO  NOISE  RATIO  FOR  CASE  1 


Signal  to  noise  ratio  for  various  values  of  cutoff 
frequency  and  anti-aliasing  filter  order  evaluated 
at  the  Nyquist  frequency  (fN)  for  a flat  input  spec 
trum. 


Signal  to  Alias  Noise  Ratios (dB) 


Filter  Order 

2 

(12dB/oct) 

4 

(24dB/oct) 

6 

(36dB/oct) 

.5 

8.4 

21.0 

33.1 

Cutoff  Frequency  (f  ) 

V 

.4 

12.3 

28.8 

44.7 

Nyquist  Frequency  ( f N) 

.315 

16.4 

36.8 

56.8 

.25 

20.5 

45.1 

69.2 

.20 

24.4 

52.9 

80.8 

.16 

28.4 

60.9 

92.9 

.125 

32.6 

69.2 

105.4 

8-6 


O 


equivalent  (as  might  be  expected)  to  adding  another  pole  to 
the  antialiasing  filter,  at  frequencies  below  the  start  of 
roll  off.  However,  above  the  frequency  where  roll  off  occurs, 
the  improvement  decreases  because  the  level  of  the  signal  is 
decreasing;  i.e., 

SN_ (n, f ) fM  2 7 

SJ5iT_?)  . tjL,  /(l+(f/fr)2)  (B.4) 

Thus  if  f <<f„  the  signal  to  noise  ratio  at  f = fXT  will 
be  the  same  as  in  the  case  of  the  flat  spectrum!  If  f is 
chosen  as  significantly  less  than  fN,  then  at  f the  signal 
to  noise  ratio  is  improved  by  20  log  (fN/f  ) . 

Thus  we  conclude  that  the  effect  of  a small  amount  of 
aliasing  is  to  introduce  a frequency-dependent  signal  to  noise 
ratio  that  depends  to  some  extent  on  the  shape  of  the  spectrum. 

For  a flat  or  falling  spectrum,  the  worst-case  signal  to  noise 
ratio  occurs  at  the  maximum  frequency  of  interest.  As  is 
expected  the  aliasing  noise  error  can  be  reduced  to  an  arbitrarily 
small  amount  by  choosing  an  anti-aliasing  filter  with  sufficiently 
steep  attenuation  characteristics,  and  setting  the  cutoff  fre- 
quency sufficiently  far  below  the  Nyquist  frequency. 
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APPENDIX  C 


EXPERIMENT  PLAN  - SUS  MINI* EXPERIMENT 

C . 1 . SITUATION 

Signal,  Underwater  Sound  (SUS)  are  widely  used  as  sources 
for  experimental  evaluation  of  long  range  propagation.  Recent 
studies  have  indicated  the  desirability  for  improving  upon 
the  accuracy  of  calibration  of  these  sources.  As  part  of  an 
analysis  of  the  feasibility  of  performing  a large  scale  experi- 
ment,. it  is  necessary  to  identify  the  inherent  variability  of 
the  source  in  terms  of  fundamental  par.ijneters. 


C-l 


C.2.  TECHNICAL  OBJECTIVES 

The  technical  objectives  of  the  program  are: 

• Identify  the  yield  and  depth  variability  inherent 
in  SUS  charge  firings  using  a statistically  signif- 
icant data  sample. 

• To  the  extent  possible,  identify  the  variability 
of  the  frequency  spectrum  of  the  surface-reflected 
wave  within  the  frequency  range  of  observation. 

• Accumulate  a data  base  from  which  additional  infor- 
mation could  be  derived  from  additional  signal  proc- 
essing. 

The  concept  of  the  experiment  that  will  satisfy  these 
objectives  is  to  set  off  a long  series  of  SUS  shots  at  the 
several  depths  of  interest,  and  record  the  resulting  acoustic 
signals  using  appropriately  positioned  and  calibrated  signal 
acquisition  systems  so  that  the  depth  of  the  shot  can  be  deter- 
mined to  acceptable  accuracy  from  the  relative  time  of  arrival 
of  direct  and  surface-reflected  signals  only.  This  data  plus 
the  bubble  pulse  frequency  of  the  shot  permits  estimation  of 
shot  yield  independent  of  the  magnitude  of  the  acoustic  signal 
received.  The  variability  of  the  surface  reflection  may  be 
identified  by  direct  comparison  of  the  spectra  of  the  direct  and 
reflected  signals  when  they  can  be  separated  in  time,  or  more 
sophisticated  signal  processing  when  they  cannot  be  separated. 
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C.  3 . TECHNICAL  PLAN 


C.3.1  General 

The  plan  by  which  data  will  be  acquired  to  fulfill  the 
technical  objectives  is  described  below. 

C-3.1.1  Overview 

The  acquisition  of  data  during  this  exercise  is  expected 
to  be  completed  within  a time  period  of  about  five  hours  plus 
time  for  deployment  and  recovery  of  data  acquisition  systems. 

The  planned  deployment  of  instrumentation  is  shown  schematically 
in  Figure  C-l,  while  a plan  view  of  the  experimental  configura- 
tion and  shot  region  are  shown  in  Figure  02.  As  shown  in 
Figure  02,  the  shots  are  fired  in  an  annular  region  of  inner 
radius  500  m and  an  outer  radius  of  1000  m.  The  signal  acquisi- 
tion systems  planned  for  use  in  this  exercise  are  the  PAR  system 
(Programmable  Acoustic  Recorder).  PAR  1,  as  deployed,  is  defined 
as  the  center  of  the  shot  circle.  It  has  deep  hydrophones 
placed  at  nominally  2750  m and  3050  m depths,  which  are  used  to 
identify  the  angle  from  vertical,  and  depth  of  the  shot.  The 
signals  will  then  be  used  to  study  the  variation  of  surface- 
reflections  at  near-vertical  incidence.  PAR  2 is  placed  at 
3000  m horizontal  range  from  PAR  1.  It  has  shallow  hydrophones 
at  500  m and  800  m depths,  which  will  be  used  as  backups  for 
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Figure  C-l  Ins trumentation  Deployment 


depth  determination  and  to  study  surface  reflection  effects 
at  near-grazing  incidence.  Both  hydrophone  arrays  will  be  pro- 
vided with  transponders  that  will  be  used  to  identify  ship 
position  relative  to  the  array,  as  well  as  relative  positions 
of  the  arrays.  The  SUS  are  dropped  in  the  annular  region  so 
as  to  produce  a nominally  constant  orientation  to  PAR  1,  while 
providing  a nominal  2/1  range  variation  to  PAR  2.  A total  of 
100  shots  are  to  be  fired  at  each  of  three  depths  (18  m,  91  m 
and  274  m respectively) . A circular  transit  of  the  source 
area  at  a nominal  speed  of  6 knots,  using  a radius  that  places 
the  ship  near  the  outer  perimeter  of  the  source  area  , will 
require  approximately  30  minutes.  The  source  area  layout  and 
experiment  timing  are  based  on  the  assumption  that  the  ship 
will  transit  the  area  in  30  minutes  (corresponding  to  a turn 
rate  of  0.2°/sec,  for  6.0  kts  ship  speed  and  a circle  radius 
of  870  m) . During  each  circle  maneuver,  33  shots  will  be 
fired,  requiring  a tota)  of  nine  circle  maneuvers  to  be  performed. 
The  same  depth  shot  will  be  used  for  all  shots  on  each  circle, 
but  depths  will  be  changed  between  circles  (i.e.,  33  18  tn,  fol- 
lowed by  33  91  m,  followed  by  33  274  m) . Use  of  approximately 
90°  of  the  annular  sector  is  not  desirable  because  the  presence 
of  the  ship  will  distort  the  surface-reflection  to  PAR  2 and 
use  of  the  transponder  will  interfere  with  data  acquisition. 
Therefore  each  of  the  33  shots  will  be  fired  at  intervals  of 
40  +5  seconds  so  that  the  total  shot  firing  time  per  circle  is 
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22  +0.5  minutes.  The  angular  sector  not  used  depends  upon 
where  the  firing  mechanism  is  located  on  the  ship  and  the  ship 
path.  The  sector  indicated  in  Figure  2 (i.e.  270°  to  360°  from 
the  line  between  PAR  1 and  2)  is  appropriate  for  launching  of 
SUS  from  the  ship  stern  at  a bearing  of  135°  to  180°  from  ship 
centerline , and  ship  circling  to  the  right.  The  last  90°  of 
the  circle  is  to  be  used  to  verify  ship  tracking  relative  to 
PAR  1,  and  time  of  crossing  the  axis  of  PAR  1 - PAR  2,  as 
determined  from  the  signals  received  from  the  transponders. 
During  this  time  interval,  the  PAR  units  will  also  be  performing 
internal  calibration  programs,  on  alternate  circles. 


C.3.1.2  Schedule 


C . 3 . 1 . 2 . 1 Ship  Schedule 
Time  from  first  SUS 

-1200  - 0000  Deck  calibration  and  deploy  PAR  1 and  PAR  2. 

PAR  in-water  pre-cal.  Determine  PAR  1 and 
PAR  2 positions,  and  array  depths.  Deploy 
deep  XBT, 

Start  first  circle.  Drop  33  19  m SUS 
Complete  first  circle,  check  position  re  PAR  1 
and  axis  crossing  re  PAR  2. 

Start  second  circle.  Drop  33  91  ® SUS 
Complete  second  circle,  check  position  re  PAR  1 
and  axis  crossing  re  PAR  2.  Deploy  XBT. 


0000  - 0023 
0023  - 0030 

0330  - 0053 
0053  - 0100 
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0057  - 0058 
0100  - 0430 
0445  - 0500 


PAR  internal  calibration  cycle 
Repeat  above  for  total  of  9 circles. 
PAR  internal  calibration. 


C. 3. 1.2. 2 PAR  Calibration  Schedule 
Time  Restive  to  first  shot 
-0030  - 0014  In-water  Pre-cal 

22  - 30  second  duration  sine  wave  cali- 
brations (see  frequency  table. 

Table  C-l)  (11  min) 

6-30  second  duration  square  wave 
calibrations,  6 Hz,  11  Hz,  17  Hz, 


21  Hz,  34  Hz,  41  Hz  (3  min) 

1-1  minute  3.33  Hz  rectangular  wave 

calibration  (1  min) 

1-1  minute  preamplifier  calibration  (1  min) 

0057  - 0058  \ 

015?  - 0158  I 1 minute,  10  Hz  square  wave  cal  at  preamps 
0257  - 0258  V to  all  hydrophones. 

0357  - 0350  1 
0457  - 0458  / 

0530  - 0546  jln-Water  Post-Cai;  repeat  above  Pre-CalJ 
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TABLE  C-l 


PAR  Frequency  Calibration  Table 
Uso  Two  Synthesizers  Sine-Sine 


Frequency  Table  — Hz 

Cal.  Number  Synthesizer  1 Synthesizer  2 


1 

14 

157 

2 

20 

163 

3 

27 

170 

4 

34 

177 

5 

41 

184 

6 

48 

191 

7 

54 

197 

8 

61 

204 

9 

68 

211 

10 

75 

218 

il 

82 

225 

12 

88 

2 31 

13 

95 

238 

14 

102 

245 

15 

109 

252 

16 

116 

259 

17 

12  3 

266 

18 

129 

272 

19 

136 

279 

20 

143 

286 

21 

150 

293 

22 

157 

300 

NOTE:  Frequency  for  synthesizer  2 is  frequency  for 

synthesizer  1 plus  143  Hz.  Frequencies  may  be 
reordered,  or  traded  to  reduce  proqranuninq  size. 
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C.3.2  Experimental  Operations 


C . 3 . 2 . 1 Signal  Measurements 

Signal  will  be  acquired  on  two  PAR  units  of  two  hydrophones 
each,  with  data  from  each  hydrophone  to  be  recorded  on  four 
(4)  tape  recorder  channels.  The  channel  gains  are  to  be  set 
with  the  anticipation  of  a nominal  peak  pressure  in  the  water 
of  199  dB  re  luPa.  The  nominal  peak  pressure  expected  in  the 
bC-300  Hz  band  is  188  dB  re  luPa.  The  four  tape  recorder 
channels  are  to  be  set  with  nominal  peak  record  sensitivities 
in  the  0 - 300  Hz  passband  region  as  follows  (dB  re  luPa) : 

a)  196  dB 

b)  190  dB 

c)  184  dB 

d)  178  dB 

Those  record  channels  should  all  be  sot  for  300  Hz  low-pass 
filter  and  phase  reference  marker  on. 

The  calibration  schedule  for  the  PAR  micro-processor  should 
be  programmed  according  to  Section  3, 1.2.2.  Calibration  peak 
signal  amplitudes  for  sine  and  square  waves  should  be  set  to 
correspond  to  2 dB  below  nominal  full  scale  on  the  tape  recorder 
(i.e.,  for  1,414  V full  scale  recorder  response  (H  harmonic 
distortion),  set  peak  sine  and  square  wave  voltages  to  1.12 
V)  . 
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C.3.2.2  Environmental  Measurements 


The  required  environmental  and  test  geometry  measurements 
include: 

a.  Depth  of  hydrophone  transponders  determined  to  within 
+3  m. 

b.  Horizontal  range  between  PAR  1 and  PAR  2 to  +15  m. 

c.  XBT  recordings  of  velocity  profile  to  depth  of  deep 

array  twice  (before  and  after  experiment)  and  four 
shallow  XBT  (hourly)  during  the  data  acquisition. 

d.  Vertical  separation  between  hydrophones  on  each  array 
shall  be  determined  to  +0.3  m. 

e.  Determination  of  the  fall  time  of  each  shot  shall 

be  determined  on  board  to  +0.01  min  or  +0.5  sec. 

f.  Ship  horizontal  range  to  PAR  1 should  be  determined 
to  a specific  radius  to  within  +10G  m.  Radius  to 
be  identified  on  board. 

g.  Charts  showing  nominal  location  of  each  shot  (bearinq 
+5°,  range  +100  m relative  to  PAR  1)  shall  be  main- 
tained. 

h.  Meteorological  conditions  and  sea  state  shall  be 
logged  hourly,  and  significant  changes  noted  when  they 


C.3.2.3  Ship  Operations 

As  noted  in  the  above  sections  the  required  ship  operations 
are  as  follows: 

• Deploy  PAR  units 

• Deploy  deep  XBT 

• Maneuver  to  identify  PAR  array  depths  and  horizontal 
separation  using  installed  transponder. 

• Perform  nine  30-minute  circling  maneuvers  around 

PAR  1,  deploying  33  SUS  shots  per  circle,  and  shallow 
XBT  on  every  alternate  circle. 


C.4  DATA  ANALYSIS 


C.4.1  Discussion 

The  objectives  of  the  SUS  mini-experiment  are  to  provide 
data  on;  a)  the  variability  of  yield  and  depth  and  b)  the  vari- 
ability of  surface-reflected  rays  from  SUS  explosions.  The 
data  analysis  deals  separately  with  these  two  objectives. 

The  experiment  is  designed  so  that  yield  and  depth  can 
be  determined  from  time  domain  measurements  only.  These  specific 
measurements  include;  a)  difference  in  time  of  arrival  of  the 
direct  signal  at  two  different  vertical  hydrophones;  and  for 
each  hydrophone,  b)  difference  in  time  of  arrival  between  direct 
and  surface-reflected  signal,  and  c)  time  between  shock  and 
bubble  pulse  peak  (bubble  pulse  period) . Data  analysis 
will  consist  of  reviewing  the  multi-rack  recordings  of  the 
signals  from  four  hydrophones  of  300  shots,  and  identifying 
the  above  time  parameters.  These  parameters  will  then  be  used 
as  inputs  to  provide  estimates  of  depth  and  yield  for  each  shot. 
Statistical  analysis  of  depth  and  yield  estimates  will  then 
provide  parameters  that  describe  shot  variability.  Representative 
spectra  of  selected  shots  will  be  compared  to  identify  the  signif- 
icance of  the  yield  and  depth  variations  on  the  resulting  source 
level. 

The  second  experiment  objective,  to  study  surface-reflection 
variability,  is  accomplished  by  comparison  of  spectra  of  isolated 
direct  and  surface-reflected  signals  from  individual  shots.  The 
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proposed  form  of  this  comparison  is  the  difference  between  energy 
spectra  of  direct  and  surface-reflected  signals.  The  average 
value  of  this  difference  over  all  similar  shots  represents  a 
measure  of  the  surface-reflection  loss.  The  standard  deviation 
(i.e.  rms  value)  of  this  difference  is  a measure  of  the  surface- 
reflection  variability.  This  approach  is  sufficient  when  the 
signals  are  separated  in  time,  but  is  not  applicable  for  shallow 
(18  m)  shots  where  significant  overlap  occurs.  In  this  case 
advanced  signal  processing  techniques  to  deconvolve  the  direct 
and  surface-reflected  signals  are  required  to  separate  them 
for  the  above  analysis.  The  principal  approach  for  attacking 
this  problem  is  expected  to  involve  correlation  analysis  to 
identify  portions  of  the  signal  that  are  unique  to  the  direct 
and  surface-reflected  components. 

C.4.2  Task  Statements 

1.  Analyze  recordings  of  300  shots  from  SUS  mini- 
experiment to  determine  yield  and  depth  variability  and 
statistical  measures  of  the  variation.  Identify  the 
impact  of  this  variation  on  the  acoustic  signal. 

2,  Analyze  recordings  to  identity  variations  in  surface- 
reflected  signals.  Such  an  analysis  is  to  include  identi- 
fication of  mean  and  variation  of  the  difference  in  direct- 
and  surface-reflected  energy  spectra. 
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MEMORANDUM  FOR  DISTRIBUTION  LIST 

Subj : DECLASSIFICATION  OF  LONG  RANGE  ACOUSTIC  PROPAGATION  PROJECT 
(LRAPP)  DOCUMENTS 
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number  of  classified  LRAPP  documents. 
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